Organize data

Merge output data from multiple batches, select metabolite columns, add participant info and other data

acg_batch1 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021-16-15.csv")

acg_batch2 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021-16-21.csv")

acg_LCM_corr <- bind_rows(acg_batch1, acg_batch2)
anyDuplicated(acg_LCM_corr$File_ID) #check for duplicated values
## [1] 106
acg_LCM_corr<-unique(acg_LCM_corr) #remove duplicate rows
acg_LCM_corr["4", "Cho_GPC_PCh_conc_molar"]<-NA #Cho conc invalid for PS10008_PS14_046_2386 due to spike in spectrum
acg_LCM_corr["4", "Cho_GPC_PCh_conc_molal"]<-NA #Cho conc invalid for PS10008_PS14_046_2386 due to spike in spectrum

write.csv(acg_LCM_corr, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021_all.csv", row.names = FALSE)

acg_molal <- acg_LCM_corr %>%
  dplyr::select(File_ID, fGM, fWM, fCSF, NAA_conc_molal, Cre_PCr_conc_molal, Cho_GPC_PCh_conc_molal, Ins_conc_molal, Glu_conc_molal, Glu_Gln_conc_molal)

acg_demo <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_analysis_subs.csv")

acg_dat <- left_join(acg_demo, acg_molal, by = "File_ID")
write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_n113.csv", row.names = FALSE)

lag_batch1 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-13-2021-15-18.csv")

lag_batch2 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021-11-41.csv")

lag_batch3 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021-12-18.csv")

lag_LCM_corr <- bind_rows(lag_batch1, lag_batch2, lag_batch3)
write.csv(lag_LCM_corr, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021_all.csv", row.names = FALSE)

lag_molal <- lag_LCM_corr %>%
  dplyr::select(File_ID, fGM, fWM, fCSF, NAA_conc_molal, Cre_PCr_conc_molal, Cho_GPC_PCh_conc_molal, Ins_conc_molal, Glu_conc_molal, Glu_Gln_conc_molal)

lag_demo <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_analysis_subs.csv")

lag_dat <- left_join(lag_demo, lag_molal, by = "File_ID")
lag_dat$subj_id<-as.character(lag_dat$subj_id) #make subj_id character type

write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_n321.csv", row.names = FALSE)

##Check data distribution & outliers Create plots to check data distribution and identify outliers

ACG

plot_age <- ggplot(acg_dat, aes(mri_age_y))
plot_age + geom_density()

#age skewed

plot_naa <- ggplot(acg_dat, aes(NAA_conc_molal))
plot_naa + geom_density()

plot_naa + geom_boxplot()

rosnerTest(acg_dat$NAA_conc_molal, k = 4)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4 
## 2.759147 2.825357 2.694361 2.721635 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 4 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 
## 3.425263 3.422302 3.419309 3.416284 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 4 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 13.36240 15.05227 14.57780 13.83339 14.71206 16.20111 15.44740 16.98446
##   [9] 15.67072 15.95324 14.87793 12.74216 16.29728 14.71133 14.66048 13.28760
##  [17] 14.59532 13.55773 16.85293 15.33699 12.45017 15.04043 15.27049 11.24177
##  [25] 14.98401 14.18774 15.83415 13.91788 17.10858 12.09629 14.70242 14.40906
##  [33] 16.46456 15.58993 15.03289 14.84971 11.30097 13.87568 15.24533 15.99486
##  [41] 14.53574 12.16795 14.20751 14.18561 16.47436 14.52629 15.00648 16.20561
##  [49] 13.22605 13.57437 14.05939 13.81423 17.33147 16.31303 16.29984 16.85100
##  [57] 15.24641 14.93994 15.14031 15.13625 18.29604 13.89531 13.42919 15.16179
##  [65] 13.87656 15.48760 12.64660 16.01636 13.64857 16.43270 15.59565 11.64927
##  [73] 15.46789 14.18371 12.72496 16.73666 16.33378 15.15593 15.86369 15.28939
##  [81] 15.31529 14.67805 13.38081 15.01442 13.01089 17.11696 14.54570 15.04389
##  [89] 14.64002 15.57552 15.21758 15.80866 15.16132 14.92979 14.18939 16.96225
##  [97] 13.86479 15.05179 15.22074 15.65731 15.51013 14.39392 15.03669 16.39137
## [105] 14.47850 14.65614 15.94476 16.78988 16.07591 13.54665 13.48426 14.73463
## [113] 15.68117
## 
## $data.name
## [1] "acg_dat$NAA_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 14.88986 1.322180 11.24177      24 2.759147   3.425263   FALSE
## 2 1 14.92243 1.281772 11.30097      37 2.825357   3.422302   FALSE
## 3 2 14.95506 1.239990 18.29604      61 2.694361   3.419309   FALSE
## 4 3 14.92469 1.203472 11.64927      72 2.721635   3.416284   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_cre <- ggplot(acg_dat, aes(Cre_PCr_conc_molal))
plot_cre + geom_density()

plot_cre + geom_boxplot()

rosnerTest(acg_dat$Cre_PCr_conc_molal, k = 5)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5 
## 3.394290 2.833607 2.707716 2.578893 2.649738 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 5 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 
## 3.425263 3.422302 3.419309 3.416284 3.413225 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 5 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 10.869437 10.445542 14.138594 11.979460 11.842544 12.567296 12.183610
##   [8] 13.315617 13.532807 13.910100 12.102060  9.609688 12.566391 11.689508
##  [15] 12.650725 10.976876 11.888348 11.823341 12.975541 12.087174 10.745566
##  [22] 12.005758 13.958950 11.012768 13.437536 12.359373 12.418194 10.775703
##  [29] 13.584172 11.251205 11.993304 12.190118 14.756858 13.780875 12.273252
##  [36] 11.243397  8.406874 11.279114 11.602557 13.761404 12.044563 10.276615
##  [43] 10.918920 10.896575 11.867174 12.189638 12.580644 12.364280 11.517322
##  [50] 10.442494 11.903800 12.566843 12.411654 12.733123 12.852863 14.066514
##  [57] 10.742053 11.977272 12.534036 12.467809 13.778825 11.508522 11.740032
##  [64] 11.569008 10.891687 12.474281 11.395342 12.716648 12.033736 13.646045
##  [71] 12.873670  9.632426 11.901341 12.610355  9.187364 13.087309 12.160354
##  [78] 11.717870 12.392428 11.411606 12.417766 11.214681 11.894967 11.689931
##  [85] 11.069854 13.648131 12.837410 11.764882 12.329981 11.666319 12.825568
##  [92] 13.032058 12.630243 12.355787 11.038718  9.731136 12.044142 11.736949
##  [99] 11.772225 11.427900 12.588592 11.074387 11.520514 12.986814 12.563423
## [106] 12.334515 13.095183 12.401466 11.143206 11.414014 11.499070 12.331493
## [113] 11.900989
## 
## $data.name
## [1] "acg_dat$Cre_PCr_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i      SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 12.03590 1.0691568  8.406874      37 3.394290   3.425263   FALSE
## 2 1 12.06830 1.0167044  9.187364      75 2.833607   3.422302   FALSE
## 3 2 12.09426 0.9833373 14.756858      33 2.707716   3.419309   FALSE
## 4 3 12.07005 0.9540396  9.609688      12 2.578893   3.416284   FALSE
## 5 4 12.09263 0.9284692  9.632426      72 2.649738   3.413225   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_cho <- ggplot(acg_dat, aes(Cho_GPC_PCh_conc_molal))
plot_cho + geom_density()
## Warning: Removed 1 rows containing non-finite values (stat_density).

plot_cho + geom_boxplot()
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

rosnerTest(acg_dat$Cho_GPC_PCh_conc_molal, k = 4)
## Warning in rosnerTest(acg_dat$Cho_GPC_PCh_conc_molal, k = 4): 1 observations
## with NA/NaN/Inf in 'x' removed.
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4 
## 2.692059 2.486294 2.542571 2.587426 
## 
## $sample.size
## [1] 112
## 
## $parameters
## k 
## 4 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 
## 3.422302 3.419309 3.416284 3.413225 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 4 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 3.054728 2.338141 2.483788 2.615586 2.904744 3.029351 2.945605 2.946224
##   [9] 2.469310 2.647417 2.329234 2.525650 2.964839 2.377024 2.316960 2.388015
##  [17] 2.587590 2.964150 2.987201 1.884924 2.294818 2.160623 1.723047 2.273437
##  [25] 2.486538 2.364455 2.110030 3.042283 2.145223 2.626855 2.657139 3.307897
##  [33] 2.993478 2.528009 2.694001 2.751861 2.169544 2.576529 2.777803 2.298135
##  [41] 2.093485 2.245856 1.887058 3.332061 2.683945 2.960541 2.306091 2.400287
##  [49] 2.262987 2.373107 2.418340 2.428692 2.651362 2.724707 2.627568 3.213338
##  [57] 2.541659 2.436647 2.503022 3.168533 2.273151 2.230455 2.575725 2.463493
##  [65] 2.964749 2.427496 2.382196 2.398624 2.592440 2.075228 1.967159 2.210616
##  [73] 2.178031 1.572637 2.265866 2.193475 2.435427 2.749229 2.738993 2.096637
##  [81] 2.376787 2.150157 2.031053 1.736687 2.549619 2.032952 2.153716 2.354125
##  [89] 2.585598 2.624180 2.597158 2.503254 2.679817 2.109408 2.755517 2.618999
##  [97] 2.254582 2.379855 2.105926 2.389925 2.621861 2.349337 2.962557 2.620008
## [105] 2.643999 2.559238 3.322683 2.403398 2.301871 2.505683 2.762249 2.704160
## 
## $data.name
## [1] "acg_dat$Cho_GPC_PCh_conc_molal"
## 
## $bad.obs
## [1] 1
## 
## $all.stats
##   i   Mean.i      SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 2.495888 0.3429536 1.572637      75 2.692059   3.422302   FALSE
## 2 1 2.504206 0.3329678 3.332061      45 2.486294   3.419309   FALSE
## 3 2 2.496680 0.3248691 3.322683     108 2.542571   3.416284   FALSE
## 4 3 2.489102 0.3164519 3.307897      33 2.587426   3.413225   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_ins <- ggplot(acg_dat, aes(Ins_conc_molal))
plot_ins + geom_density()

shapiro.test(acg_dat$Ins_conc_molal) #not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  acg_dat$Ins_conc_molal
## W = 0.96211, p-value = 0.002754
plot_ins + geom_boxplot()

rosnerTest(acg_dat$Ins_conc_molal, k = 3) #outlier value 1.471244, observation 84 (PS18_013) - spectrum showed issue, exclude from analysis 
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3 
## 3.753416 3.044570 3.051178 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 3 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 
## 3.425263 3.422302 3.419309 
## 
## $n.outliers
## [1] 1
## 
## $alternative
## [1] "Up to 3 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  7.633881  8.178435  7.218240  7.696923  8.696692  8.600324  9.402207
##   [8]  8.875539  6.779586  5.935284 10.304907  5.028532  7.058793  9.067489
##  [15]  6.700483  7.361395  8.124760  8.682546  6.553380  6.320046  7.315123
##  [22]  7.446418  7.281237  7.102747  9.505161  7.839894  7.887706  7.616099
##  [29]  7.601724  6.347168  7.912922  5.656884  9.656654  9.672370  8.341248
##  [36]  8.816154  5.196891  6.304340  6.789199  8.018694  7.594973  7.416233
##  [43]  7.406756  3.198494  2.985922  9.898020  8.079535  8.703354  7.157223
##  [50]  6.213808  8.399767  9.174211  9.656691  9.335033 10.080614  4.813291
##  [57]  9.728369  7.521052  7.441993  8.051283 10.887902  8.695335  5.872553
##  [64]  6.418561  5.242212  8.749670  7.765345  8.185832  7.950928  8.394652
##  [71]  9.583041  6.202420  7.829477  6.209769  4.285017  8.677590  7.775045
##  [78]  9.334317  8.684385  7.236647  7.742767  6.219638  9.105159  1.471244
##  [85]  5.218826  9.407651  7.634559  6.064850  8.785683  3.770173  7.618745
##  [92]  8.763466  7.123140  6.399259  6.804541  5.855502  9.364210  6.203701
##  [99] 10.565682  9.856594  7.453487  7.728481  7.846516  9.947994  8.640216
## [106]  7.758280  8.476271  8.672477  6.437713  5.739577  6.673856  7.596796
## [113]  7.229168
## 
## $data.name
## [1] "acg_dat$Ins_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 7.588864 1.629881 1.471244      84 3.753416   3.425263    TRUE
## 2 1 7.643485 1.529793 2.985922      45 3.044570   3.422302   FALSE
## 3 2 7.685445 1.470564 3.198494      44 3.051178   3.419309   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
acg_dat["84", "Ins_conc_molal"] <- NA #replace outlier value with NA to exclude from analysis

plot_glx <- ggplot(acg_dat, aes(Glu_Gln_conc_molal))
plot_glx + geom_density()

plot_glx + geom_boxplot()

rosnerTest(acg_dat$Glu_Gln_conc_molal, k = 1)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1 
## 2.840029 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 1 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 
## 3.425263 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 1 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 22.16208 24.54895 30.79902 25.15445 23.12144 23.45996 27.95219 28.85494
##   [9] 28.47293 30.83024 24.02170 20.83895 26.81225 25.81936 25.83210 28.15351
##  [17] 28.69984 24.45281 25.87441 26.10568 24.90431 27.69072 29.16664 25.73952
##  [25] 31.00567 25.29233 29.35298 28.81394 30.45098 26.42377 25.25458 30.67838
##  [33] 29.65833 34.25604 26.24353 28.06087 21.93668 25.00839 27.05607 29.92424
##  [41] 27.12366 23.65348 26.82576 24.86136 26.92752 23.98651 26.21693 28.06806
##  [49] 26.26637 27.44238 23.98526 25.40089 22.59948 23.46421 28.23542 33.70878
##  [57] 26.99835 25.59893 26.95525 25.77541 32.34599 24.66776 25.46001 25.71687
##  [65] 29.73266 25.54233 27.95401 29.18637 26.09421 27.74434 27.70571 21.88957
##  [73] 26.72328 30.44184 21.32831 31.44932 28.73199 25.08766 29.57015 24.10930
##  [81] 27.84161 24.85586 29.21009 20.93551 23.34771 27.38610 26.35670 26.10330
##  [89] 28.90132 26.67329 24.94565 29.44010 29.32746 26.23233 24.32566 24.39736
##  [97] 26.88103 30.91320 28.56651 27.47953 25.89977 24.22166 25.89081 28.12079
## [105] 29.77257 28.04112 29.21096 29.61668 27.08694 28.10033 24.05989 26.66364
## [113] 26.45077
## 
## $data.name
## [1] "acg_dat$Glu_Gln_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 26.82935 2.615006 34.25604      34 2.840029   3.425263   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113.csv", row.names = FALSE)

LAG

plot_age <- ggplot(lag_dat, aes(mri_age_y))
plot_age + geom_density()

#age skewed

plot_naa <- ggplot(lag_dat, aes(NAA_conc_molal))
plot_naa + geom_density()

plot_naa + geom_boxplot()

rosnerTest(lag_dat$NAA_conc_molal, k = 10) #outliers val. 19.207, obs. 296 & Val. 10.148, obs. 306 - some issues in spectra, exclude from analysis for NAA
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 4.467031 4.044357 3.599759 3.613937 3.039885 2.720506 2.647000 2.651977 
##      R.9     R.10 
## 2.660497 2.662447 
## 
## $sample.size
## [1] 321
## 
## $parameters
##  k 
## 10 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
##  lambda.1  lambda.2  lambda.3  lambda.4  lambda.5  lambda.6  lambda.7  lambda.8 
##  3.742579  3.741706  3.740829  3.739949  3.739067  3.738181  3.737291  3.736399 
##  lambda.9 lambda.10 
##  3.735503  3.734604 
## 
## $n.outliers
## [1] 2
## 
## $alternative
## [1] "Up to 10 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 14.16473 13.06076 14.35969 13.27563 16.92063 13.92045 14.90883 14.96705
##   [9] 13.87538 13.79928 15.43119 13.74363 13.42774 14.32941 15.29816 15.36852
##  [17] 14.96343 14.01462 15.11779 16.29892 15.28430 15.17137 15.34158 15.28575
##  [25] 13.61378 16.76759 13.63702 14.41263 13.69061 15.64950 14.60942 15.25717
##  [33] 15.68290 14.24180 15.96102 13.19084 14.22355 14.45688 13.89590 16.68967
##  [41] 15.14115 13.88211 14.56337 14.77687 14.00542 14.12098 13.36743 13.70315
##  [49] 14.59050 15.55276 13.33791 14.13365 12.08635 14.08027 13.90090 13.63249
##  [57] 15.72428 11.86871 14.17563 11.96658 12.37485 14.57472 13.54554 14.35760
##  [65] 13.86325 16.89273 14.68644 16.67785 15.93255 16.85334 14.93659 14.21829
##  [73] 13.99673 14.23697 13.86957 13.86406 14.36716 14.41372 12.45632 14.76312
##  [81] 14.75832 14.92673 13.88514 14.24089 15.42366 14.86462 14.32717 14.96667
##  [89] 13.72426 14.16948 15.16894 13.93828 12.92563 13.94977 13.64160 14.55102
##  [97] 14.51660 12.52392 14.21054 13.93670 16.19979 14.45536 14.10612 13.35969
## [105] 16.58862 14.41649 14.80811 14.26277 12.37256 15.44893 13.43728 13.57859
## [113] 13.72929 14.86513 11.95501 14.39891 14.11323 14.14168 14.73091 14.40429
## [121] 13.19741 13.71443 14.26528 14.92876 14.67063 14.35362 15.12298 13.95432
## [129] 13.68703 13.44733 14.01869 14.75568 14.74986 11.38458 12.08162 11.74260
## [137] 13.93366 14.16181 14.38636 14.40623 14.15278 14.17590 15.46643 14.09616
## [145] 13.62893 14.14878 13.54441 13.09883 12.79401 14.64345 14.02714 14.59131
## [153] 14.51512 15.41640 13.41648 13.54853 13.78080 13.79989 13.51420 14.77851
## [161] 15.60295 13.13390 13.99763 14.27699 15.00866 15.01634 12.09852 14.98452
## [169] 13.60870 16.34165 11.91384 12.33377 14.09992 15.06924 13.87378 13.47019
## [177] 14.02295 18.06089 14.25601 14.37746 13.24330 15.14502 13.95585 13.43714
## [185] 14.30145 16.30323 13.44700 14.70962 13.78765 13.43206 14.86922 14.73623
## [193] 14.54106 13.30455 14.61897 13.91250 15.44818 13.00756 14.21959 16.58017
## [201] 14.30707 15.14338 12.83493 13.96649 15.41109 13.77502 12.24289 15.19071
## [209] 13.32940 13.11120 15.00722 15.15742 13.96143 14.70099 14.27321 15.14983
## [217] 14.79365 15.19557 14.84745 13.19832 15.12751 13.95088 13.28456 14.90679
## [225] 12.94820 13.28282 14.39416 13.76430 14.91352 15.55618 15.20787 14.89958
## [233] 14.31622 14.53646 15.29391 13.51980 15.05685 14.73951 15.01094 15.61748
## [241] 14.09588 14.63147 15.59956 15.33848 14.07436 14.30913 13.94217 16.47833
## [249] 15.29907 13.85813 15.73137 15.08295 16.15169 14.28221 15.30824 14.47596
## [257] 14.47137 14.28477 15.33016 14.75400 14.54703 14.66072 13.49926 14.47973
## [265] 13.83620 14.83537 15.05792 14.18100 13.68113 13.28288 14.53224 13.92267
## [273] 14.80805 13.25470 14.60733 13.77756 15.36693 16.13917 14.36568 14.54767
## [281] 14.51644 16.31083 15.95003 13.41067 15.50889 13.79439 14.53462 14.83200
## [289] 14.42306 13.93431 13.79584 14.91428 12.66542 14.98849 15.38162 19.20747
## [297] 13.74790 17.99338 14.95991 13.93438 15.42154 13.27326 13.14415 14.47594
## [305] 13.99067 10.14812 13.73230 14.91913 15.85376 13.23125 13.40927 14.90803
## [313] 14.30314 14.54287 14.33013 16.55184 14.94880 14.30460 15.31501 14.46473
## [321] 16.31675
## 
## $data.name
## [1] "lag_dat$NAA_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##    i   Mean.i      SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1  0 14.39149 1.0781156 19.20747     296 4.467031   3.742579    TRUE
## 2  1 14.37644 1.0454867 10.14812     306 4.044357   3.741706    TRUE
## 3  2 14.38970 1.0198436 18.06089     178 3.599759   3.740829   FALSE
## 4  3 14.37815 1.0003560 17.99338     298 3.613937   3.739949   FALSE
## 5  4 14.36675 0.9810139 11.38458     134 3.039885   3.739067   FALSE
## 6  5 14.37619 0.9680502 11.74260     136 2.720506   3.738181   FALSE
## 7  6 14.38455 0.9580956 16.92063       5 2.647000   3.737291   FALSE
## 8  7 14.37647 0.9488234 16.89273      66 2.651977   3.736399   FALSE
## 9  8 14.36843 0.9395707 11.86871      58 2.660497   3.735503   FALSE
## 10 9 14.37644 0.9303092 16.85334      70 2.662447   3.734604   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["296", "NAA_conc_molal"]<- NA
lag_dat["306", "NAA_conc_molal"]<- NA

plot_cre <- ggplot(lag_dat, aes(Cre_PCr_conc_molal))
plot_cre + geom_density()

plot_cre + geom_boxplot()

rosnerTest(lag_dat$Cre_PCr_conc_molal, k = 10) #outliers val. 14.829, obs. 241 (PS17_047) - structure in resid. & val. 14.09, obs. 178 (PS17_045, spectrum looked fine), excluding only PS17_047 for now
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 4.835902 4.220654 3.332561 3.278883 3.286432 3.307998 2.911341 2.895521 
##      R.9     R.10 
## 2.767415 2.775535 
## 
## $sample.size
## [1] 321
## 
## $parameters
##  k 
## 10 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
##  lambda.1  lambda.2  lambda.3  lambda.4  lambda.5  lambda.6  lambda.7  lambda.8 
##  3.742579  3.741706  3.740829  3.739949  3.739067  3.738181  3.737291  3.736399 
##  lambda.9 lambda.10 
##  3.735503  3.734604 
## 
## $n.outliers
## [1] 2
## 
## $alternative
## [1] "Up to 10 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  9.979902  8.752796  9.874363  8.961297 11.510180 10.121568 10.418593
##   [8] 11.028922 10.171413  9.929791 10.871830 11.064228 10.111309  9.319128
##  [15] 10.662947 10.890644 10.197139  9.919287 10.426231 11.325982 10.557574
##  [22] 11.088082 11.039646 11.172931  8.936929 12.109674 10.000368 10.420059
##  [29]  9.572731 10.511121  9.984628 10.695378 11.039338  9.822119 11.097416
##  [36]  9.239262  9.742404 10.192047 10.845538 12.117826 10.389433  9.535353
##  [43] 10.010295 10.823639  9.876324 10.019161  8.608303  9.442680 11.444015
##  [50] 11.419644  9.408142 10.139106  8.483950  9.792171 10.032519  8.968350
##  [57] 11.562417  9.156420 10.285179  7.449145 10.159106 10.623228 10.245682
##  [64] 11.514351 10.182653 11.786704 11.560836 12.034078 11.166365 11.463697
##  [71] 10.397742 10.176308 11.796639 10.733740 10.174742 10.465407  9.824296
##  [78] 11.420357  7.832529 10.284534 12.513826 11.576531 10.541650 10.808215
##  [85] 10.483803 10.041962 10.382300  9.976776 10.765842 10.520082 11.422196
##  [92] 10.843472  9.968723 10.781936  9.841352  9.930816 11.494602  9.908886
##  [99] 10.765624 10.122441 11.914466 10.403873 10.056177  9.179081 11.252406
## [106] 10.011672 10.666794  9.981113  9.225683 10.367639 10.879950 11.335061
## [113]  9.425501 10.725975 10.241112  9.317008  9.889557  9.560765 10.593153
## [120] 11.263291  8.935236 10.137273 10.048855 11.211108 10.481657 10.028783
## [127] 11.460008  9.844337 10.207849  8.876822  9.557376  9.418484  9.631575
## [134] 11.004550  7.414051  8.984416  9.798072 10.236338 10.150536 10.108983
## [141] 10.653763 10.137920 10.541536  9.525275  9.886573  9.896683  9.185657
## [148]  9.642602  7.281217  9.968306 10.715037 10.739283 11.039381 11.250002
## [155]  9.759444  9.643949  9.675426  9.621214  9.141950 10.356422 10.788507
## [162]  9.362209  9.822128  9.645301 10.157630 10.381588  8.913720 11.109833
## [169]  9.890247 12.488886  9.115647  9.814255 10.627685 11.734993 10.295594
## [176] 10.173568  9.988731 14.090723 10.829617 10.247883 10.533434 10.999808
## [183] 10.001382  9.229915 10.354994 11.296388  9.994843 10.904378  9.345731
## [190]  9.618365  9.759703 10.530182 10.844392  8.895843 10.110129 11.423088
## [197] 11.096328  9.958072  9.931007 10.921273  9.154313 10.482202 10.423407
## [204] 10.629254 12.212044  9.913344  8.305599 10.446897  9.841379  9.379922
## [211]  9.844787 10.968766 10.801476  9.891839 10.632212  9.671280 10.872873
## [218] 10.666877 10.645958  9.493899 13.103099  9.512697 10.276565  9.846960
## [225]  9.773848  9.929345 11.034653 10.095542 10.422810 10.640722 10.407355
## [232] 11.450386 10.284795 10.179735 11.263605 10.183927 10.257883 10.937157
## [239] 10.529431 11.012918 14.829071  9.133491 11.531961 10.173887  9.610385
## [246]  9.349643  9.324653 10.926724 11.044735  8.859819 11.228731 10.816240
## [253] 11.271292  9.861827 10.852278 10.764773  9.380198  9.481786 10.781017
## [260]  8.985104  9.888393  9.434081  8.870117  9.122209  9.542512  9.889146
## [267] 10.769356  9.616033  8.966342  9.255494 10.764599  9.897261  9.427714
## [274]  8.163350 10.017130  9.265501 10.944724 12.004282  9.764853  9.856451
## [281] 10.020483 12.009948 10.044769  9.667192 11.197137  9.964094 10.035066
## [288] 10.812708  9.728669  9.918824  8.843802 10.202045  9.183452  9.822267
## [295] 10.790509 12.454846  8.641534  9.507122 10.361671  9.646740 11.011202
## [302]  9.075410 10.166484  9.737467  9.615290  7.882414 10.377305 10.224366
## [309] 12.449600 10.083783  9.260131 11.009901  9.914846  9.659798 10.610323
## [316] 11.561746 11.077058  9.291690 10.160005  9.934510 11.934339
## 
## $data.name
## [1] "lag_dat$Cre_PCr_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##    i   Mean.i      SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1  0 10.26105 0.9446048 14.829071     241 4.835902   3.742579    TRUE
## 2  1 10.24678 0.9107457 14.090723     178 4.220654   3.741706    TRUE
## 3  2 10.23473 0.8862592  7.281217     149 3.332561   3.740829   FALSE
## 4  3 10.24402 0.8719682 13.103099     221 3.278883   3.739949   FALSE
## 5  4 10.23500 0.8583617  7.414051     135 3.286432   3.739067   FALSE
## 6  5 10.24393 0.8448557  7.449145      60 3.307998   3.738181   FALSE
## 7  6 10.25280 0.8313244  7.832529      79 2.911341   3.737291   FALSE
## 8  7 10.26051 0.8213002  7.882414     306 2.895521   3.736399   FALSE
## 9  8 10.26810 0.8114875 12.513826      81 2.767415   3.735503   FALSE
## 10 9 10.26091 0.8027211 12.488886     170 2.775535   3.734604   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["241", "Cre_PCr_conc_molal"]<-NA  #replace outlier value with NA to exclude from analysis

plot_cho <- ggplot(lag_dat, aes(Cho_GPC_PCh_conc_molal))
plot_cho + geom_density()

plot_cho + geom_boxplot()

rosnerTest(lag_dat$Cho_GPC_PCh_conc_molal, k = 8)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 3.694967 3.121637 3.061660 3.062209 2.979166 3.001536 2.871646 2.733239 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 8 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 8 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 2.068181 2.062467 2.386391 2.355035 2.498191 2.256110 2.077815 2.425292
##   [9] 2.398770 1.895982 1.865260 2.184158 2.004899 2.128662 2.236830 2.378977
##  [17] 2.355458 2.266506 2.404117 2.307009 2.052566 2.364941 2.267604 1.940308
##  [25] 2.587651 2.716570 2.469831 2.547177 2.087747 2.170207 2.232363 2.241430
##  [33] 2.441839 2.264953 2.333662 2.359801 2.270399 2.241648 1.920730 2.467142
##  [41] 2.242001 2.227010 2.229581 2.053798 1.939191 2.364334 2.678152 2.362113
##  [49] 1.676867 2.166236 2.306276 2.211050 2.178445 2.361600 2.344510 2.268087
##  [57] 2.667157 2.059386 2.034573 2.045022 2.214431 2.995506 1.846508 2.183946
##  [65] 2.072392 2.429630 2.330074 2.491408 2.265040 2.154392 2.411001 2.046897
##  [73] 2.534412 2.027012 2.376839 2.222266 2.306644 2.292419 2.108096 2.172672
##  [81] 1.839589 1.859195 1.443755 1.640698 2.420001 2.337525 2.238790 2.240345
##  [89] 2.290795 2.399979 2.350376 2.151479 2.142910 2.110996 1.985190 1.928934
##  [97] 2.304258 1.765061 2.137595 1.985591 2.781876 2.123475 1.756957 1.700318
## [105] 1.975316 1.974756 1.999963 2.031005 1.687589 1.708882 2.296279 2.182991
## [113] 2.276755 2.101586 1.837847 2.371711 2.065515 2.499424 2.363551 2.223515
## [121] 2.232809 2.045991 1.947255 2.339091 2.193248 2.528316 2.691920 2.777161
## [129] 2.451074 2.763753 2.377191 2.160428 2.215103 2.214312 2.533792 2.384939
## [137] 2.179386 2.552211 2.673232 2.476304 1.625876 1.926042 2.202001 2.204096
## [145] 2.344069 2.152159 1.895155 1.710908 1.401447 1.946252 2.047116 2.456831
## [153] 2.086612 2.201666 2.190181 2.294616 2.581968 2.225911 2.498999 2.000618
## [161] 2.163682 1.960245 2.204252 2.159811 2.097381 2.254865 1.844443 1.598729
## [169] 1.951488 3.208685 1.836756 1.950559 1.992121 2.281328 2.008449 2.599790
## [177] 2.317712 3.033330 2.698335 2.332814 2.583764 2.234143 2.191218 2.264957
## [185] 1.977277 2.067888 2.356329 2.431505 2.779803 2.150897 2.566469 2.101205
## [193] 2.159059 2.164490 2.152499 2.489624 2.791537 1.673365 2.193421 2.761471
## [201] 1.959300 2.022674 1.970196 2.366537 2.774857 2.129388 1.836375 2.360163
## [209] 2.476086 2.665161 2.258604 2.298965 2.012883 2.201029 2.158646 2.353116
## [217] 1.946114 1.912304 1.692442 2.307008 2.449170 2.132806 2.045527 2.155850
## [225] 2.166616 2.440956 1.783126 2.272129 2.595674 2.284913 2.155251 2.400964
## [233] 2.048543 2.045160 2.356748 2.222524 2.346513 2.378440 2.281038 1.906127
## [241] 1.725583 2.479767 2.262261 2.082256 2.047599 2.138599 2.252565 2.262486
## [249] 1.881557 2.322610 1.489680 2.225178 2.118903 2.390375 2.437536 2.400945
## [257] 1.832693 2.290250 2.007843 2.346441 2.089229 2.428740 2.391081 2.732220
## [265] 1.971298 2.235255 2.116196 1.941129 2.009743 1.999019 1.927678 1.994412
## [273] 2.372282 2.491822 1.886355 2.109659 2.028203 1.807961 2.439615 2.152157
## [281] 1.892138 2.163111 2.376896 1.914835 2.734225 1.944843 1.756721 2.135482
## [289] 2.101386 1.768588 1.974516 1.975407 2.167757 2.112266 2.510332 2.521732
## [297] 2.431779 2.959724 2.046249 1.986009 2.300110 2.368720 2.151772 2.070465
## [305] 2.457353 2.243696 2.105342 2.108232 2.180180 2.349512 1.911659 2.150353
## [313] 2.011103 2.243254 2.417149 1.800195 2.190437 2.875345 2.116343 2.319437
## [321] 2.422916
## 
## $data.name
## [1] "lag_dat$Cho_GPC_PCh_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i      SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 2.208213 0.2707661 3.208685     170 3.694967   3.742579   FALSE
## 2 1 2.205087 0.2653235 3.033330     178 3.121637   3.741706   FALSE
## 3 2 2.202490 0.2616371 1.401447     149 3.061660   3.740829   FALSE
## 4 3 2.205009 0.2581460 2.995506      62 3.062209   3.739949   FALSE
## 5 4 2.202516 0.2546890 1.443755      83 2.979166   3.739067   FALSE
## 6 5 2.204917 0.2514736 2.959724     298 3.001536   3.738181   FALSE
## 7 6 2.202521 0.2482341 1.489680     251 2.871646   3.737291   FALSE
## 8 7 2.204791 0.2453333 2.875345     318 2.733239   3.736399   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_ins <- ggplot(lag_dat, aes(Ins_conc_molal))
plot_ins + geom_density()

plot_ins + geom_boxplot()

rosnerTest(lag_dat$Ins_conc_molal, k = 8) #outliers: val. 14.077, obs. 241 (PS17_047) - structure in resid; val. 11.777, obs. 221 (PS16_014) - spectrum fine, keep for now; val. 1.542, obs. 197 (PS18_028) - noisy spectrum, spiky residuals, exclude; val. 2.04, obs. 62 (PS14_111) - spectrum fine, keep for now
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 5.981441 4.463363 4.077869 3.761599 3.384729 3.303295 3.148162 3.105675 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 8 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399 
## 
## $n.outliers
## [1] 4
## 
## $alternative
## [1] "Up to 8 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  8.084411  7.069186  6.913031  6.915494  7.726417  8.046174  6.937152
##   [8]  4.793633  5.952276  6.984233  6.400936  4.413732  7.176021  7.448535
##  [15]  7.547442  6.541306  8.400250  5.957678  6.965664  7.688786  6.245160
##  [22]  7.360921  7.925104  4.710055  6.035369  7.937845  7.166536  7.269577
##  [29]  5.736145  7.751891  6.217281  7.121735  7.138640  6.104017  7.736566
##  [36]  6.539681  4.958196  7.275083  6.485777  6.598573  7.153802  7.352086
##  [43]  7.119035  6.146843  6.038556  6.638197  4.629802  6.187281  6.553798
##  [50]  5.945866  7.474194  5.559932  5.743412  6.163849  8.347824  4.432658
##  [57]  8.169688  5.045280  8.873237  6.213139  6.400722  2.040275  8.158702
##  [64]  6.760986  6.468052  7.273556  6.742605  7.552488  7.843304  5.876394
##  [71]  8.162408  4.963944  7.945325  5.440451  4.854219  4.834565  5.839066
##  [78]  5.839637  5.890621  5.871859  6.963936  5.579543  5.190653  7.062975
##  [85]  6.065434  7.046003  6.865074  6.210505  6.024185  6.121049  8.133551
##  [92]  6.863480  6.537634  7.457231  5.008241  6.892333  5.936366  3.826364
##  [99]  7.075620  5.567096  5.248858  5.491554  4.529573  7.011630  5.896171
## [106]  6.268198  5.292489  5.663161  6.753982  5.643647  6.788746  5.437125
## [113]  3.078293  7.511207  2.568354  8.286250  5.309958  6.053735  8.092627
## [120]  7.122038  6.509101  6.564424  5.780997  6.770400  5.727577  7.829891
## [127]  6.206182  7.258719  6.558572  4.787960  6.132452  6.029208  6.663446
## [134]  7.630692  4.608510  5.731585  5.435869  5.209567  2.734199  5.857696
## [141]  5.638088  6.135153  6.888234  5.236533  5.655180  4.948126  4.992038
## [148]  6.025793  2.972647  6.870427  7.728083  7.800580  7.334580  6.676194
## [155]  6.170230  5.196197  7.073851  6.127912  7.232895  4.960469  6.612268
## [162]  6.088183  5.798318  6.543746  6.997835  6.987522  4.250557  4.637685
## [169]  6.778603  7.502499  4.679992  5.270903  6.367636  6.301188  7.918660
## [176]  6.533953  7.497948  8.801777  6.797100  5.264895  5.014240  6.141243
## [183]  5.472028  3.885218  4.816043  6.248540  4.626087  6.007857  6.538200
## [190]  7.586651  7.411192  7.320786  7.682231  7.120925  6.686811  7.373729
## [197]  1.541994  6.314608  6.840831  7.095049  6.617076  6.056406  5.882616
## [204]  4.775951  5.731427  5.979244  6.418407  6.735628  5.243289  5.869558
## [211]  5.226573  6.275003  7.093997  7.541310  5.731832  4.649731  6.150065
## [218]  6.428658  7.426624  6.729085 11.777066  5.591387  3.080430  8.027587
## [225]  7.066068  6.661947  6.600256  6.955873  6.949031  6.824263  7.605754
## [232]  6.414457  6.094616  6.805794  5.865261  8.270934  4.935241  6.818562
## [239]  5.441384  7.674150 14.076805  4.939637  6.137708  5.300423  6.306576
## [246]  6.245176  7.663984  7.103391  8.335725  3.649002  9.145931  7.122485
## [253]  8.331023  7.199990  8.067478  7.557525  5.521439  6.036198  6.082461
## [260]  7.524328  5.926927  5.705055  5.399996  6.232379  5.707127  4.605042
## [267]  4.333050  5.607343  5.940738  5.740379  6.177392  6.685727  5.210419
## [274]  6.865329  6.749236  5.464544  6.496419  6.705573  7.274507  6.221882
## [281]  5.994352  6.059500  7.660503  6.660801  7.439851  6.777698  5.282780
## [288]  8.773125  6.173129  5.673978  5.639679  5.174471  8.615918  5.774767
## [295]  5.411721  7.832813  7.925781  4.996084  6.276000  6.161567  6.484734
## [302]  5.594141  6.793109  6.401047  9.077790  8.562826  5.674753  6.368745
## [309]  8.425749  6.658481  4.433254  5.482190  7.149805  6.079902  6.139307
## [316]  7.111624  5.971685  6.120964  6.571466  5.418973  6.379001
## 
## $data.name
## [1] "lag_dat$Ins_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 6.383077 1.286267 14.076805     241 5.981441   3.742579    TRUE
## 2 1 6.359034 1.213890 11.777066     221 4.463363   3.741706    TRUE
## 3 2 6.342050 1.177099  1.541994     197 4.077869   3.740829    TRUE
## 4 3 6.357144 1.147615  2.040275      62 3.761599   3.739949    TRUE
## 5 4 6.370762 1.123401  2.568354     115 3.384729   3.739067   FALSE
## 6 5 6.382795 1.104533  2.734199     139 3.303295   3.738181   FALSE
## 7 6 6.394378 1.086898  2.972647     149 3.148162   3.737291   FALSE
## 8 7 6.405275 1.071259  3.078293     113 3.105675   3.736399   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["241", "Ins_conc_molal"] <- NA #replace outlier value with NA to exclude from analysis
lag_dat["197", "Ins_conc_molal"] <- NA

plot_glx <- ggplot(lag_dat, aes(Glu_Gln_conc_molal))
plot_glx + geom_density()

plot_glx + geom_boxplot()

rosnerTest(lag_dat$Glu_Gln_conc_molal, k = 7)#outliers: val. 10.09, obs. 241 (PS17_047) - structure in resid; val. 31.2899, obs. 296 (PS21_030) - big spike in resid. around 2.1, exclude.
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7 
## 4.954047 4.247538 3.396313 3.439741 3.463175 3.267721 3.232088 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 7 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 
## 
## $n.outliers
## [1] 2
## 
## $alternative
## [1] "Up to 7 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 23.58311 20.72833 21.43465 20.79838 24.55542 19.38008 22.13533 21.98148
##   [9] 18.58332 21.80094 23.39061 22.39400 22.25994 19.84839 22.62240 19.74196
##  [17] 19.55718 20.18644 20.93614 24.01329 22.14911 23.96545 23.49753 23.68093
##  [25] 20.83139 25.55075 18.83744 23.02848 20.10324 24.84118 23.11824 25.62522
##  [33] 25.14532 20.60338 23.10183 18.68727 19.90887 22.63726 23.44172 23.97287
##  [41] 20.79684 21.03860 21.41871 22.38334 19.07457 19.80466 21.54623 21.54133
##  [49] 22.80842 25.04145 20.65555 22.02043 18.52301 21.98565 22.88652 20.99602
##  [57] 23.64194 20.45706 24.44658 22.40116 19.22839 24.77315 20.40702 22.32472
##  [65] 22.79321 23.27919 23.66026 23.59194 23.28915 23.15967 20.72270 22.14974
##  [73] 23.87191 27.09158 23.12072 25.11318 22.08527 21.17324 14.25503 21.44083
##  [81] 23.08444 21.47161 23.79592 26.39922 24.52352 23.18206 22.57951 22.53444
##  [89] 24.24539 23.38651 21.39081 23.50524 21.84539 23.16617 20.66627 20.68771
##  [97] 21.33913 23.82111 24.17549 20.43359 24.17777 21.80255 21.79162 21.34246
## [105] 26.43216 21.02534 21.61161 22.57479 21.58113 22.47585 21.89183 22.86451
## [113] 21.63056 23.83333 23.16833 22.80824 21.74170 19.79860 22.68433 25.43932
## [121] 20.46600 19.55287 21.71921 23.91575 19.94617 22.54020 25.12373 20.26601
## [129] 20.42848 18.40115 22.31740 21.34145 21.18639 24.46723 17.55526 19.17975
## [137] 22.14220 21.22080 20.62138 19.29726 22.33080 22.36491 23.64798 24.04583
## [145] 25.02751 22.02316 20.84404 22.82850 23.75262 19.93399 21.93649 21.46368
## [153] 23.51970 23.79151 22.29261 22.55887 19.88319 23.52534 18.26039 22.44651
## [161] 23.87196 21.33171 21.19221 21.21896 25.13735 21.72890 21.00404 24.38403
## [169] 22.77448 24.07955 20.92792 22.79597 21.53526 25.58587 20.23221 20.38536
## [177] 19.64011 29.12643 21.85437 22.81675 24.27197 20.09271 21.28381 21.20096
## [185] 20.61092 26.55674 20.24095 22.22341 21.17598 21.05662 20.70603 22.17497
## [193] 22.22738 17.77073 21.62541 23.02400 23.16615 17.38515 20.64300 21.97048
## [201] 23.41849 22.42124 20.67899 21.35918 23.92414 18.75710 19.72431 22.35249
## [209] 23.13971 19.66645 21.97753 21.86935 22.15373 18.04793 20.63749 19.81250
## [217] 22.91987 21.99328 21.86022 19.23943 22.18618 21.97414 22.10059 21.96999
## [225] 21.38803 21.08812 21.12026 20.43893 20.01495 22.48505 23.44854 23.27656
## [233] 22.72471 21.03119 22.44160 19.23947 24.24070 21.72010 21.19354 21.24039
## [241] 10.09497 19.32810 23.68178 22.11865 18.91616 19.59606 19.02554 22.92979
## [249] 26.99382 18.91880 22.36158 21.71252 22.51005 20.90044 23.69281 22.60270
## [257] 20.69658 19.00693 19.97580 15.08181 17.70438 17.08767 17.40514 16.78328
## [265] 20.40818 19.57305 21.57031 21.73287 18.10301 16.64465 21.45075 20.54401
## [273] 19.32920 14.36555 21.66594 17.20659 21.57909 24.57946 17.49328 18.35123
## [281] 21.22156 24.11680 22.51873 19.19853 23.08553 19.03798 23.38476 21.59376
## [289] 19.09657 21.28636 18.16122 19.70762 18.41625 18.73419 22.65583 31.28994
## [297] 18.89292 17.32671 18.83142 19.60466 24.66522 18.35557 22.92964 18.49429
## [305] 21.56271 18.98307 20.17174 24.36413 25.12361 21.76663 19.89145 21.01574
## [313] 20.06260 25.41609 26.15708 23.69646 23.29001 20.54746 23.65148 20.44966
## [321] 28.51150
## 
## $data.name
## [1] "lag_dat$Glu_Gln_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 21.68857 2.340229 10.09497     241 4.954047   3.742579    TRUE
## 2 1 21.72480 2.251926 31.28994     296 4.247538   3.741706    TRUE
## 3 2 21.69482 2.190549 14.25503      79 3.396313   3.740829   FALSE
## 4 3 21.71821 2.153715 29.12643     178 3.439741   3.739949   FALSE
## 5 4 21.69484 2.116351 14.36555     274 3.463175   3.739067   FALSE
## 6 5 21.71804 2.078963 28.51150     321 3.267721   3.738181   FALSE
## 7 6 21.69647 2.046559 15.08181     260 3.232088   3.737291   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["241", "Glu_Gln_conc_molal"] <- NA
lag_dat["296", "Glu_Gln_conc_molal"] <- NA

#exclude PS17_047 from all analyses due to outlier status on multiple metabolites and issues in spectrum (noisy, structure in resid)
lag_dat["241", "NAA_conc_molal"] <- NA
lag_dat["241", "Cre_PCr_conc_molal"] <- NA
lag_dat["241", "Cho_GPC_PCh_conc_molal"] <- NA

write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320.csv", row.names = FALSE)

##Descriptive stats

#ACG - count individual subject N and N of females
acg_sub_sex <- acg_dat %>%
  dplyr::select(subj_id, female)
acg_unique_subs<-unique(acg_sub_sex)
length(acg_unique_subs$subj_id)
## [1] 102
sum(acg_unique_subs$female)
## [1] 47
describe(acg_dat)
##                        vars   n    mean      sd  median trimmed     mad     min
## study_code*               1 113   57.00   32.76   57.00   57.00   41.51    1.00
## subj_id*                  2 113   51.42   28.89   54.00   51.46   35.58    1.00
## spec_id                   3 113 5072.42 3034.48 3244.00 4708.09 1254.28 2058.00
## date_scan*                4 113   49.21   28.71   50.00   49.32   37.06    1.00
## hand*                     5 113    4.54    1.04    5.00    4.78    0.00    1.00
## female                    6 113    0.47    0.50    0.00    0.46    0.00    0.00
## mri_age_y                 7 113    4.06    1.10    3.68    3.97    1.00    2.34
## folder*                   8 113   57.00   32.76   57.00   57.00   41.51    1.00
## File_ID*                  9 113   57.00   32.76   57.00   57.00   41.51    1.00
## fGM                      10 113    0.77    0.05    0.78    0.78    0.05    0.55
## fWM                      11 113    0.04    0.03    0.03    0.03    0.02    0.00
## fCSF                     12 113    0.19    0.05    0.18    0.19    0.04    0.09
## NAA_conc_molal           13 113   14.89    1.32   15.03   14.94    1.25   11.24
## Cre_PCr_conc_molal       14 113   12.04    1.07   12.04   12.06    0.84    8.41
## Cho_GPC_PCh_conc_molal   15 112    2.50    0.34    2.48    2.49    0.30    1.57
## Ins_conc_molal           16 112    7.64    1.53    7.71    7.72    1.46    2.99
## Glu_conc_molal           17 113   22.38    2.25   22.50   22.46    2.01   16.03
## Glu_Gln_conc_molal       18 113   26.83    2.62   26.72   26.81    2.64   20.84
##                             max   range  skew kurtosis     se
## study_code*              113.00  112.00  0.00    -1.23   3.08
## subj_id*                 102.00  101.00 -0.05    -1.15   2.72
## spec_id                11448.00 9390.00  0.87    -0.81 285.46
## date_scan*                97.00   96.00 -0.04    -1.29   2.70
## hand*                      6.00    5.00 -2.07     3.38   0.10
## female                     1.00    1.00  0.12    -2.00   0.05
## mri_age_y                  7.36    5.03  0.72    -0.49   0.10
## folder*                  113.00  112.00  0.00    -1.23   3.08
## File_ID*                 113.00  112.00  0.00    -1.23   3.08
## fGM                        0.87    0.32 -1.14     1.59   0.01
## fWM                        0.29    0.29  4.07    25.28   0.00
## fCSF                       0.36    0.27  0.90     0.92   0.00
## NAA_conc_molal            18.30    7.05 -0.37     0.21   0.12
## Cre_PCr_conc_molal        14.76    6.35 -0.36     0.83   0.10
## Cho_GPC_PCh_conc_molal     3.33    1.76  0.16     0.08   0.03
## Ins_conc_molal            10.89    7.90 -0.51     0.35   0.14
## Glu_conc_molal            28.20   12.17 -0.31     0.40   0.21
## Glu_Gln_conc_molal        34.26   13.42  0.13    -0.01   0.25
#LAG - count individual subject N and N of females
lag_sub_sex <- lag_dat %>%
  dplyr::select(subj_id, female)
lag_unique_subs<-unique(lag_sub_sex)
length(lag_unique_subs$subj_id)
## [1] 95
sum(lag_unique_subs$female)
## [1] 47
describe(lag_dat)
##                        vars   n    mean      sd  median trimmed     mad     min
## study_code*               1 321  160.94   92.72  161.00  161.00  118.61    1.00
## subj_id*                  2 321   40.27   24.55   41.00   39.16   28.17    1.00
## spec_id                   3 320 7744.26 3536.93 6379.50 7565.04 3846.61 3118.00
## date_scan*                4 321  143.73   82.33  143.00  143.92  103.78    1.00
## hand*                     5 321    4.64    0.90    5.00    4.88    0.00    1.00
## female                    6 321    0.47    0.50    0.00    0.46    0.00    0.00
## mri_age_y                 7 321    5.95    1.90    5.39    5.79    1.79    2.41
## folder*                   8 321  161.00   92.81  161.00  161.00  118.61    1.00
## File_ID*                  9 321  161.00   92.81  161.00  161.00  118.61    1.00
## fGM                      10 321    0.62    0.10    0.64    0.63    0.10    0.30
## fWM                      11 321    0.35    0.11    0.34    0.34    0.11    0.00
## fCSF                     12 321    0.03    0.02    0.03    0.03    0.01    0.00
## NAA_conc_molal           13 318   14.39    1.02   14.36   14.37    0.88   11.38
## Cre_PCr_conc_molal       14 320   10.25    0.91   10.17   10.23    0.83    7.28
## Cho_GPC_PCh_conc_molal   15 320    2.21    0.27    2.20    2.20    0.24    1.40
## Ins_conc_molal           16 319    6.37    1.19    6.37    6.39    1.08    2.04
## Glu_conc_molal           17 321   17.88    1.90   17.87   17.88    1.73   10.63
## Glu_Gln_conc_molal       18 319   21.69    2.19   21.74   21.72    2.02   14.26
##                             max    range  skew kurtosis     se
## study_code*              320.00   319.00  0.00    -1.21   5.17
## subj_id*                  95.00    94.00  0.28    -0.83   1.37
## spec_id                14205.00 11087.00  0.35    -1.37 197.72
## date_scan*               284.00   283.00  0.00    -1.19   4.60
## hand*                      6.00     5.00 -2.35     4.48   0.05
## female                     1.00     1.00  0.13    -1.99   0.03
## mri_age_y                 11.13     8.72  0.66    -0.46   0.11
## folder*                  321.00   320.00  0.00    -1.21   5.18
## File_ID*                 321.00   320.00  0.00    -1.21   5.18
## fGM                        0.80     0.50 -0.75     0.40   0.01
## fWM                        0.66     0.65  0.43     0.20   0.01
## fCSF                       0.21     0.21  3.12    16.87   0.00
## NAA_conc_molal            18.06     6.68  0.20     0.95   0.06
## Cre_PCr_conc_molal        14.09     6.81  0.17     1.33   0.05
## Cho_GPC_PCh_conc_molal     3.21     1.81  0.22     0.80   0.02
## Ins_conc_molal            11.78     9.74 -0.13     1.69   0.07
## Glu_conc_molal            24.87    14.24 -0.07     1.22   0.11
## Glu_Gln_conc_molal        29.13    14.87 -0.12     0.74   0.12
#Calculate total Ns from both datasets combined
all_dat<-bind_rows(acg_dat, lag_dat)
length(unique(all_dat$File_ID))
## [1] 434
length(unique(all_dat$subj_id))
## [1] 127
all_sub_sex <- all_dat %>%
  dplyr::select(subj_id, female)
all_unique_subs<-unique(all_sub_sex)
sum(all_unique_subs$female)
## [1] 62
describe(all_dat)
##                        vars   n    mean      sd  median trimmed     mad     min
## study_code*               1 434  216.94  124.93  216.50  216.97  160.12    1.00
## subj_id*                  2 434   53.87   33.15   54.00   52.20   38.55    1.00
## spec_id                   3 433 7046.99 3606.23 5832.00 6818.88 3964.47 2058.00
## date_scan*                4 434  184.27  105.60  183.00  184.45  133.43    1.00
## hand*                     5 434    4.62    0.93    5.00    4.86    0.00    1.00
## female                    6 434    0.47    0.50    0.00    0.46    0.00    0.00
## mri_age_y                 7 434    5.46    1.92    4.95    5.28    1.68    2.34
## folder*                   8 434  217.04  125.00  217.50  217.05  160.12    1.00
## File_ID*                  9 434  217.50  125.43  217.50  217.50  160.86    1.00
## fGM                      10 434    0.66    0.11    0.67    0.67    0.11    0.30
## fWM                      11 434    0.27    0.17    0.29    0.27    0.16    0.00
## fCSF                     12 434    0.07    0.08    0.03    0.06    0.03    0.00
## NAA_conc_molal           13 431   14.52    1.13   14.48   14.51    0.99   11.24
## Cre_PCr_conc_molal       14 433   10.71    1.24   10.54   10.65    1.10    7.28
## Cho_GPC_PCh_conc_molal   15 432    2.28    0.32    2.26    2.27    0.26    1.40
## Ins_conc_molal           16 431    6.70    1.40    6.66    6.68    1.28    2.04
## Glu_conc_molal           17 434   19.05    2.81   18.53   18.89    2.53   10.63
## Glu_Gln_conc_molal       18 432   23.04    3.23   22.54   22.85    2.77   14.26
##                             max    range  skew kurtosis     se
## study_code*              432.00   431.00  0.00    -1.21   6.00
## subj_id*                 127.00   126.00  0.33    -0.78   1.59
## spec_id                14205.00 12147.00  0.45    -1.19 173.30
## date_scan*               366.00   365.00  0.00    -1.18   5.07
## hand*                      6.00     5.00 -2.29     4.27   0.04
## female                     1.00     1.00  0.13    -1.99   0.02
## mri_age_y                 11.13     8.79  0.79    -0.09   0.09
## folder*                  433.00   432.00  0.00    -1.20   6.00
## File_ID*                 434.00   433.00  0.00    -1.21   6.02
## fGM                        0.87     0.57 -0.60     0.10   0.01
## fWM                        0.66     0.66 -0.09    -0.87   0.01
## fCSF                       0.36     0.36  1.37     0.83   0.00
## NAA_conc_molal            18.30     7.05  0.10     0.57   0.05
## Cre_PCr_conc_molal        14.76     7.48  0.44     0.20   0.06
## Cho_GPC_PCh_conc_molal     3.33     1.93  0.45     0.70   0.02
## Ins_conc_molal            11.78     9.74  0.08     0.70   0.07
## Glu_conc_molal            28.20    17.57  0.49     0.13   0.13
## Glu_Gln_conc_molal        34.26    20.00  0.54     0.34   0.16

LME modelling ACG

run mixed-effects models for ACG with participant as random effect, age, sex and age x sex interaction as fixed effects NAA run on acg_data_molal_outliers_removed_n113.csv

acg_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    383.3    394.2   -187.6    375.3      109 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.53495 -0.39238  0.05273  0.38044  1.68117 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 1.1187   1.0577  
##  Residual             0.6121   0.7824  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  13.8379     0.4588 108.1511   30.16   <2e-16 ***
## mri_age_y     0.2513     0.1079 105.1446    2.33   0.0217 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.959
acg_naa_age_sex<- lmer(NAA_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    385.2    398.8   -187.6    375.2      108 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.54531 -0.40334  0.05382  0.37656  1.67013 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 1.1189   1.0578  
##  Residual             0.6113   0.7819  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  13.85437    0.46506 109.75797  29.791   <2e-16 ***
## mri_age_y     0.25366    0.10837 104.64089   2.341   0.0211 *  
## female       -0.05648    0.26075  99.22276  -0.217   0.8289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.926       
## female    -0.165 -0.098
anova(acg_naa_age,acg_naa_age_sex)
## Data: acg_dat
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_age_sex: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_naa_age        4 383.25 394.16 -187.62   375.25                     
## acg_naa_age_sex    5 385.20 398.84 -187.60   375.20 0.0469  1     0.8285
acg_naa_ageXsex<- lmer(NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    387.1    403.4   -187.5    375.1      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55254 -0.39337  0.03275  0.37793  1.69016 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 1.1179   1.0573  
##  Residual             0.6105   0.7813  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       13.6886     0.6657 106.4622  20.562   <2e-16 ***
## mri_age_y          0.2954     0.1616 107.8036   1.827   0.0704 .  
## female             0.2522     0.9254 110.1751   0.273   0.7857    
## mri_age_y:female  -0.0757     0.2178 109.0452  -0.348   0.7288    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.964              
## female      -0.719  0.694       
## mri_g_y:fml  0.716 -0.742 -0.960
anova(acg_naa_age,acg_naa_ageXsex)
## Data: acg_dat
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_ageXsex: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_naa_age        4 383.25 394.16 -187.62   375.25                     
## acg_naa_ageXsex    6 387.08 403.45 -187.54   375.08 0.1677  2     0.9196
#including sex in model does not significantly improve fit 

Cre

acg_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    340.9    351.8   -166.5    332.9      109 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12075 -0.39494  0.02556  0.37255  1.62021 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.6867   0.8287  
##  Residual             0.4877   0.6983  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  11.67366    0.38192 110.29328  30.566   <2e-16 ***
## mri_age_y     0.09221    0.08992 107.78610   1.026    0.307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.960
acg_cre_age_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cre_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    342.9    356.5   -166.4    332.9      108 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1065 -0.3932  0.0080  0.3605  1.6020 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.6866   0.8286  
##  Residual             0.4871   0.6979  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  11.65746    0.38678 111.48623  30.139   <2e-16 ***
## mri_age_y     0.08988    0.09034 107.34558   0.995    0.322    
## female        0.05586    0.21445  96.67410   0.260    0.795    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.927       
## female    -0.160 -0.100
anova(acg_cre_age,acg_cre_age_sex)
## Data: acg_dat
## Models:
## acg_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cre_age_sex: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cre_age        4 340.94 351.85 -166.47   332.94                     
## acg_cre_age_sex    5 342.87 356.51 -166.44   332.87 0.0678  1     0.7945
acg_cre_ageXsex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    343.6    360.0   -165.8    331.6      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.11740 -0.37592  0.04868  0.38312  1.59263 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.7009   0.8372  
##  Residual             0.4638   0.6810  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       12.09093    0.54699 105.10134  22.104   <2e-16 ***
## mri_age_y         -0.01931    0.13288 106.80701  -0.145    0.885    
## female            -0.77294    0.76525 111.06817  -1.010    0.315    
## mri_age_y:female   0.20325    0.18016 110.02199   1.128    0.262    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.965              
## female      -0.715  0.689       
## mri_g_y:fml  0.711 -0.738 -0.960
anova(acg_cre_age,acg_cre_ageXsex)
## Data: acg_dat
## Models:
## acg_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cre_ageXsex: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cre_age        4 340.94 351.85 -166.47   332.94                     
## acg_cre_ageXsex    6 343.61 359.98 -165.81   331.61 1.3273  2      0.515

Cho

acg_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     78.1     89.0    -35.0     70.1      108 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1325 -0.5979 -0.1286  0.6106  2.4400 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.03204  0.1790  
##  Residual             0.07887  0.2808  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.74592    0.12031 111.92574   22.82   <2e-16 ***
## mri_age_y    -0.06423    0.02842 111.94035   -2.26   0.0258 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.963
acg_cho_age_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cho_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     78.2     91.8    -34.1     68.2      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17836 -0.63945 -0.09041  0.54644  2.50032 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.03214  0.1793  
##  Residual             0.07694  0.2774  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.77074    0.12063 111.75493  22.970   <2e-16 ***
## mri_age_y    -0.06004    0.02834 111.86929  -2.118   0.0364 *  
## female       -0.09017    0.06493 103.29750  -1.389   0.1679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.931       
## female    -0.148 -0.107
anova(acg_cho_age,acg_cho_age_sex)
## Data: acg_dat
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_age_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cho_age        4 78.077 88.951 -35.039   70.077                     
## acg_cho_age_sex    5 78.166 91.759 -34.083   68.166 1.9111  1     0.1668
acg_cho_ageXsex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cho_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     79.7     96.1    -33.9     67.7      106 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1308 -0.6109 -0.1102  0.5606  2.5638 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.03297  0.1816  
##  Residual             0.07579  0.2753  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        2.69546    0.16726 107.68134  16.116   <2e-16 ***
## mri_age_y         -0.04107    0.04072 109.53005  -1.008    0.315    
## female             0.05961    0.23991 111.98464   0.248    0.804    
## mri_age_y:female  -0.03674    0.05661 111.93901  -0.649    0.518    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.965              
## female      -0.697  0.673       
## mri_g_y:fml  0.694 -0.719 -0.963
anova(acg_cho_age,acg_cho_ageXsex)
## Data: acg_dat
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_ageXsex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cho_age        4 78.077 88.951 -35.039   70.077                     
## acg_cho_ageXsex    6 79.747 96.058 -33.874   67.747 2.3297  2      0.312
#including sex in model does not significantly improve fit 

Ins

acg_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    416.5    427.3   -204.2    408.5      108 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.51194 -0.47558  0.09217  0.51196  1.57868 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.9268   0.9627  
##  Residual             1.3749   1.1726  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   8.4742     0.5464 111.9652  15.509   <2e-16 ***
## mri_age_y    -0.2049     0.1294 111.4483  -1.583    0.116    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.962
acg_ins_age_sex<- lmer(Ins_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    417.1    430.7   -203.5    407.1      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.47103 -0.48748  0.05729  0.54401  1.51171 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.8891   0.9429  
##  Residual             1.3819   1.1755  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   8.5595     0.5485 111.9724  15.605   <2e-16 ***
## mri_age_y    -0.1859     0.1296 111.3196  -1.435    0.154    
## female       -0.3489     0.2981  99.5526  -1.171    0.245    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.930       
## female    -0.140 -0.118
anova(acg_ins_age,acg_ins_age_sex)
## Data: acg_dat
## Models:
## acg_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_ins_age_sex: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_ins_age        4 416.45 427.33 -204.23   408.45                     
## acg_ins_age_sex    5 417.09 430.69 -203.55   407.09 1.3598  1     0.2436
acg_ins_ageXsex<- lmer(Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    418.2    434.5   -203.1    406.2      106 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.55603 -0.53421  0.07378  0.57588  1.65390 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.7615   0.8727  
##  Residual             1.4790   1.2162  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)        9.0615     0.7667 105.7180  11.819   <2e-16 ***
## mri_age_y         -0.3138     0.1882 108.1469  -1.667   0.0983 .  
## female            -1.3693     1.0920 111.9981  -1.254   0.2125    
## mri_age_y:female   0.2511     0.2588 111.8723   0.970   0.3340    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.965              
## female      -0.702  0.678       
## mri_g_y:fml  0.702 -0.727 -0.963
anova(acg_ins_age,acg_ins_ageXsex)
## Data: acg_dat
## Models:
## acg_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_ins_ageXsex: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_ins_age        4 416.45 427.33 -204.23   408.45                     
## acg_ins_ageXsex    6 418.20 434.51 -203.10   406.20 2.2549  2     0.3239

Glx

acg_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    541.1    552.0   -266.5    533.1      109 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.48535 -0.37243 -0.00651  0.43846  1.68169 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.706    2.169   
##  Residual             2.330    1.526   
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  25.7502     0.9206 106.5717  27.972   <2e-16 ***
## mri_age_y     0.2612     0.2163 103.0785   1.207     0.23    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.959
acg_glx_age_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    543.0    556.6   -266.5    533.0      108 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.47879 -0.39442 -0.02857  0.42851  1.70104 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.692    2.166   
##  Residual             2.335    1.528   
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  25.7002     0.9334 108.5734  27.535   <2e-16 ***
## mri_age_y     0.2542     0.2173 102.5984   1.170    0.245    
## female        0.1700     0.5257  98.8828   0.323    0.747    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.925       
## female    -0.167 -0.097
anova(acg_glx_age,acg_glx_age_sex)
## Data: acg_dat
## Models:
## acg_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_glx_age_sex: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_glx_age        4 541.09 552.00 -266.55   533.09                     
## acg_glx_age_sex    5 542.99 556.62 -266.49   532.99 0.1045  1     0.7465
acg_glx_ageXsex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    542.5    558.9   -265.2    530.5      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.53333 -0.43763  0.07324  0.38166  1.63409 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.573    2.139   
##  Residual             2.297    1.515   
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       27.2150     1.3265 106.5825  20.516   <2e-16 ***
## mri_age_y         -0.1271     0.3220 107.8777  -0.395    0.694    
## female            -2.6328     1.8384 109.4802  -1.432    0.155    
## mri_age_y:female   0.6875     0.4325 108.2254   1.590    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.964              
## female      -0.722  0.696       
## mri_g_y:fml  0.718 -0.745 -0.959
anova(acg_glx_age,acg_glx_ageXsex)
## Data: acg_dat
## Models:
## acg_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_glx_ageXsex: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## acg_glx_age        4 541.09 552.00 -266.55   533.09                    
## acg_glx_ageXsex    6 542.49 558.85 -265.25   530.49 2.603  2     0.2721

Plot results ACG

Create scatter plot with regression line based on model predictions NAA

acg_dat <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113.csv")

#pull intercept and slope from model for regression line
fixef.acg_naa<-fixef(acg_naa_age)

#plot
plot_acg_naa_age<-ggplot(acg_dat, aes(x=mri_age_y, y=NAA_conc_molal,color=as.factor(female)))
plot_acg_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.acg_naa[[1]], slope= fixef.acg_naa[[2]])) + 
  theme_classic()

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_naa_age_viridis.png")
## Saving 7 x 5 in image

Cho

#pull intercept and slope from model for regression line
fixef.acg_cho<-fixef(acg_cho_age)

#plot
plot_acg_cho_age<-ggplot(acg_dat, aes(x=mri_age_y, y=Cho_GPC_PCh_conc_molal,color=as.factor(female)))
plot_acg_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "tCho (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.acg_cho[[1]], slope= fixef.acg_cho[[2]])) + 
  theme_classic()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_cho_age_viridis.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 row(s) containing missing values (geom_path).

Ins

#pull intercept and slope from model for regression line
fixef.acg_ins<-fixef(acg_ins_age)

#plot
plot_acg_ins_age<-ggplot(acg_dat, aes(x=mri_age_y, y=Ins_conc_molal,color=as.factor(female)))
plot_acg_ins_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "Ins (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.acg_ins[[1]], slope= fixef.acg_ins[[2]])) + 
  theme_classic()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_ins_age_viridis.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 row(s) containing missing values (geom_path).

LME modelling LAG

run mixed-effects models for lag with participant as random effect, age, sex and age x sex interaction as fixed effects NAA, substantial longitudinal data so we include random slope run on lag_data_molal_outliers_removed_n320.csv

#basic model effect of age on intercept
lag_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##      896      911     -444      888      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4313 -0.5792 -0.0185  0.5440  3.6263 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1641   0.4051  
##  Residual             0.8272   0.9095  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  13.82375    0.18588 309.72511  74.369  < 2e-16 ***
## mri_age_y     0.09699    0.02902 316.51634   3.343  0.00093 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.929
lag_naa_age_sex<- lmer(NAA_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    898.0    916.8   -444.0    888.0      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4286 -0.5791 -0.0185  0.5452  3.6231 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1637   0.4046  
##  Residual             0.8275   0.9097  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  13.827819   0.194854 290.918032  70.965  < 2e-16 ***
## mri_age_y     0.097086   0.029038 316.589446   3.343 0.000927 ***
## female       -0.009878   0.137557  93.706753  -0.072 0.942903    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.874       
## female    -0.300 -0.039
anova(lag_naa_age,lag_naa_age_sex)
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age_sex: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_naa_age        4 895.96 911.01 -443.98   887.96                     
## lag_naa_age_sex    5 897.95 916.76 -443.98   887.95 0.0051  1      0.943
lag_naa_ageXsex<- lmer(NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    899.9    922.5   -444.0    887.9      312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4372 -0.5859 -0.0225  0.5510  3.6185 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1645   0.4056  
##  Residual             0.8268   0.9093  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       13.79473    0.24812 314.11810  55.598  < 2e-16 ***
## mri_age_y          0.10272    0.03910 313.83715   2.627  0.00903 ** 
## female             0.06577    0.37459 308.18552   0.176  0.86075    
## mri_age_y:female  -0.01265    0.05838 316.98288  -0.217  0.82864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.924              
## female      -0.662  0.612       
## mri_g_y:fml  0.619 -0.670 -0.930
anova(lag_naa_age,lag_naa_ageXsex)
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_ageXsex: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_naa_age        4 895.96 911.01 -443.98   887.96                     
## lag_naa_ageXsex    6 899.91 922.48 -443.95   887.91 0.0518  2     0.9744
#including sex in model does not significantly improve fit 

#test quadratic effect of age
#create age^2  variable
lag_dat$mri_age_y2 <- (lag_dat$mri_age_y)^2

lag_naa_age2<- lmer(NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age2) #age^2 p=.0447
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    893.9    912.7   -442.0    883.9      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5451 -0.5819 -0.0270  0.4941  3.5330 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1700   0.4123  
##  Residual             0.8119   0.9011  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  12.70181    0.58588 310.95949  21.680   <2e-16 ***
## mri_age_y     0.47130    0.18749 303.10334   2.514   0.0125 *  
## mri_age_y2   -0.02836    0.01402 297.51463  -2.023   0.0440 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) mr_g_y
## mri_age_y  -0.983       
## mri_age_y2  0.949 -0.988
anova(lag_naa_age, lag_naa_age2) #fit improved p=.044
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age2: NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## lag_naa_age     4 895.96 911.01 -443.98   887.96                       
## lag_naa_age2    5 893.91 912.72 -441.96   883.91 4.0455  1    0.04429 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cre

lag_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    848.2    863.3   -420.1    840.2      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9764 -0.5307 -0.1261  0.5525  4.1548 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08477  0.2912  
##  Residual             0.73683  0.8584  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.49912    0.16805 308.49608  62.476   <2e-16 ***
## mri_age_y    -0.04183    0.02646 319.87112  -1.581    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.937
lag_cre_age_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    849.7    868.6   -419.9    839.7      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9528 -0.5658 -0.1208  0.5434  4.1925 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08478  0.2912  
##  Residual             0.73557  0.8577  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.46482    0.17478 289.38001  59.873   <2e-16 ***
## mri_age_y    -0.04267    0.02646 319.88902  -1.612    0.108    
## female        0.08318    0.11739  84.73092   0.709    0.481    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.887       
## female    -0.277 -0.044
anova(lag_cre_age,lag_cre_age_sex)
## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_age_sex: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_cre_age        4 848.23 863.31 -420.12   840.23                     
## lag_cre_age_sex    5 849.73 868.57 -419.87   839.73 0.5017  1     0.4788
lag_cre_ageXsex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    850.0    872.6   -419.0    838.0      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9579 -0.5653 -0.0910  0.5718  4.2323 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08177  0.2860  
##  Residual             0.73339  0.8564  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       10.65304    0.22484 314.11103  47.381   <2e-16 ***
## mri_age_y         -0.07476    0.03585 318.90184  -2.086   0.0378 *  
## female            -0.33761    0.33688 307.33953  -1.002   0.3171    
## mri_age_y:female   0.07047    0.05297 319.97938   1.330   0.1844    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.935              
## female      -0.667  0.624       
## mri_g_y:fml  0.632 -0.677 -0.938
anova(lag_cre_age,lag_cre_ageXsex)
## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_ageXsex: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## lag_cre_age        4 848.23 863.31 -420.12   840.23                    
## lag_cre_ageXsex    6 849.97 872.58 -418.98   837.97 2.263  2     0.3226

Cho

lag_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     51.3     66.4    -21.7     43.3      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6539 -0.5810 -0.0562  0.4659  4.1202 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.0149   0.1221  
##  Residual             0.0561   0.2368  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.346670   0.049038 310.144986  47.854  < 2e-16 ***
## mri_age_y    -0.023315   0.007581 315.773363  -3.075  0.00229 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.921
lag_cho_age_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     51.6     70.5    -20.8     41.6      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6317 -0.5893 -0.0587  0.4688  4.1717 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01469  0.1212  
##  Residual             0.05586  0.2364  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.36786    0.05147 283.09492  46.006  < 2e-16 ***
## mri_age_y    -0.02286    0.00757 315.96555  -3.019  0.00274 ** 
## female       -0.05024    0.03808  80.73098  -1.319  0.19083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.860       
## female    -0.312 -0.046
anova(lag_cho_age,lag_cho_age_sex)
## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_age_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## lag_cho_age        4 51.345 66.419 -21.673   43.345                    
## lag_cho_age_sex    5 51.610 70.452 -20.805   41.610 1.735  1     0.1878
lag_cho_ageXsex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     49.9     72.5    -18.9     37.9      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6639 -0.5411 -0.0413  0.4642  4.0526 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01452  0.1205  
##  Residual             0.05521  0.2350  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        2.28998    0.06505 314.68723  35.206   <2e-16 ***
## mri_age_y         -0.00954    0.01019 312.12139  -0.936   0.3497    
## female             0.12488    0.09791 309.31059   1.275   0.2031    
## mri_age_y:female  -0.02931    0.01511 316.65469  -1.939   0.0533 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.916              
## female      -0.664  0.609       
## mri_g_y:fml  0.617 -0.674 -0.922
anova(lag_cho_age,lag_cho_ageXsex)
## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_ageXsex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## lag_cho_age        4 51.345 66.419 -21.673   43.345                       
## lag_cho_ageXsex    6 49.871 72.481 -18.936   37.871 5.4741  2    0.06476 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#including sex in model does not significantly improve fit 

Ins

lag_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1013.3   1028.3   -502.6   1005.3      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3984 -0.5594 -0.0111  0.5911  4.6516 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1768   0.4205  
##  Residual             1.2227   1.1057  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.28473    0.21938 307.28604  28.648   <2e-16 ***
## mri_age_y     0.01632    0.03444 318.45543   0.474    0.636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.934
lag_ins_age_sex<- lmer(Ins_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1015.3   1034.1   -502.6   1005.3      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4003 -0.5588 -0.0120  0.5901  4.6506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1765   0.4201  
##  Residual             1.2229   1.1058  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   6.286355   0.228788 284.644902  27.477   <2e-16 ***
## mri_age_y     0.016356   0.034473 318.479871   0.474    0.635    
## female       -0.003854   0.157311  78.609245  -0.024    0.981    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.881       
## female    -0.284 -0.046
anova(lag_ins_age,lag_ins_age_sex)
## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_age_sex: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## lag_ins_age        4 1013.3 1028.3 -502.64   1005.3                    
## lag_ins_age_sex    5 1015.3 1034.1 -502.64   1005.3 6e-04  1     0.9807
lag_ins_ageXsex<- lmer(Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1017.1   1039.7   -502.6   1005.1      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4248 -0.5418 -0.0121  0.5899  4.6348 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1767   0.4204  
##  Residual             1.2221   1.1055  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        6.361751   0.294417 312.731852  21.608   <2e-16 ***
## mri_age_y          0.003458   0.046822 316.620308   0.074    0.941    
## female            -0.171702   0.441468 305.644875  -0.389    0.698    
## mri_age_y:female   0.028161   0.069170 318.668940   0.407    0.684    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.930              
## female      -0.667  0.620       
## mri_g_y:fml  0.630 -0.677 -0.934
anova(lag_ins_age,lag_ins_ageXsex)
## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_ageXsex: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_ins_age        4 1013.3 1028.3 -502.64   1005.3                     
## lag_ins_ageXsex    6 1017.1 1039.7 -502.56   1005.1 0.1663  2     0.9202

Glx

lag_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
isSingular(lag_glx_age) # TRUE: model is almost/near singular (parameters are on the boundary of the feasible parameter space: random effects variance = 0)
## [1] TRUE
summary(lag_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1396.1   1411.1   -694.0   1388.1      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6700 -0.6848  0.0127  0.6729  3.5338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             4.542    2.131   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  23.23676    0.39326 319.00000  59.088  < 2e-16 ***
## mri_age_y    -0.25961    0.06309 319.00000  -4.115 4.93e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
lag_glx_age_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1396.7   1415.5   -693.4   1386.7      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7491 -0.6907  0.0192  0.7055  3.6028 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             4.523    2.127   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  23.12126    0.40478 319.00000  57.120  < 2e-16 ***
## mri_age_y    -0.26202    0.06299 319.00000  -4.160  4.1e-05 ***
## female        0.27789    0.23880 319.00000   1.164    0.245    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.915       
## female    -0.245 -0.033
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(lag_glx_age,lag_glx_age_sex)
## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_age_sex: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_glx_age        4 1396.1 1411.1 -694.04   1388.1                     
## lag_glx_age_sex    5 1396.7 1415.5 -693.36   1386.7 1.3513  1      0.245
lag_glx_ageXsex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1397.8   1420.4   -692.9   1385.8      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7059 -0.6942 -0.0013  0.6785  3.6197 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.00     0.000   
##  Residual             4.51     2.124   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       23.44862    0.52588 319.00000  44.589  < 2e-16 ***
## mri_age_y         -0.31768    0.08502 319.00000  -3.737 0.000221 ***
## female            -0.45362    0.78868 319.00000  -0.575 0.565583    
## mri_age_y:female   0.12295    0.12636 319.00000   0.973 0.331265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.951              
## female      -0.667  0.634       
## mri_g_y:fml  0.640 -0.673 -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(lag_glx_age,lag_glx_ageXsex)
## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_ageXsex: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_glx_age        4 1396.1 1411.1 -694.04   1388.1                     
## lag_glx_ageXsex    6 1397.8 1420.4 -692.89   1385.8 2.2967  2     0.3172

Plot results LAG

Create scatter plot with regression line based on model predictions NAA

lag_dat<-read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320.csv")

#pull intercept and slope from model for regression line
fixef.lag_naa<-fixef(lag_naa_age)

#plot
plot_lag_naa_age<-ggplot(lag_dat, aes(x=mri_age_y, y=NAA_conc_molal,color=as.factor(female)))
plot_lag_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_naa[[1]], slope= fixef.lag_naa[[2]])) + 
  theme_classic()
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_naa_age_viridis.png")
## Saving 7 x 5 in image
## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 3 row(s) containing missing values (geom_path).

Cho

#pull intercept and slope from model for regression line
fixef.lag_cho<-fixef(lag_cho_age)

#plot
plot_lag_cho_age<-ggplot(lag_dat, aes(x=mri_age_y, y=Cho_GPC_PCh_conc_molal,color=as.factor(female)))
plot_lag_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tCho (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_cho[[1]], slope= fixef.lag_cho[[2]])) + 
  theme_classic()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_cho_age_viridis.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 row(s) containing missing values (geom_path).

Glx

#pull intercept and slope from model for regression line
fixef.lag_glx<-fixef(lag_glx_age)

#plot
plot_lag_glx_age<-ggplot(lag_dat, aes(x=mri_age_y, y=Glu_Gln_conc_molal,color=as.factor(female)))
plot_lag_glx_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "Glx (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_glx[[1]], slope= fixef.lag_glx[[2]])) + 
  theme_classic()
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_glx_age_viridis.png")
## Saving 7 x 5 in image
## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 row(s) containing missing values (geom_path).

##Combine results plots into 1 figure Simplify plots to legend & region don’t repeat in each panel

library(patchwork)
p1_acg_naa_age<-plot_acg_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.acg_naa[[1]], slope= fixef.acg_naa[[2]])) + 
  theme_classic() +
  theme(plot.title= element_text(hjust = 0.5))

p4_lag_naa_age<-plot_lag_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  theme_classic() +
  #theme(legend.position = c(0.7, 0.95), legend.background = element_rect(fill='transparent')) +
  geom_abline(aes(intercept = fixef.lag_naa[[1]], slope= fixef.lag_naa[[2]])) +
  theme(plot.title= element_text(hjust = 0.5))


p2_acg_cho_age<-plot_acg_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "tCho (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.acg_cho[[1]], slope= fixef.acg_cho[[2]])) + 
  theme_classic() 
  

p5_lag_cho_age<-plot_lag_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "tCho (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_cho[[1]], slope= fixef.lag_cho[[2]])) + 
  theme_classic()

p3_acg_ins_age<-plot_acg_ins_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "Ins (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.acg_ins[[1]], slope= fixef.acg_ins[[2]])) + 
  theme_classic()

p6_lag_glx_age<-plot_lag_glx_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "Glx (molal)") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_glx[[1]], slope= fixef.lag_glx[[2]])) + 
  theme_classic()

results_fig<-p1_acg_naa_age + p4_lag_naa_age + p2_acg_cho_age + p5_lag_cho_age + p3_acg_ins_age + p6_lag_glx_age + plot_layout(nrow=3, guides ='collect')

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/results_fig_grid_viridis.png", results_fig, width = 2000, units = "px")
## Saving 2000 x 1500 px image
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).

##Correct for multiple comparisons Use FDR method (Benjamini & Hochberg) to correct for mutliple comparisons (5 metabolites x 2 voxels = 10 models)

#list models (voxel/metabolite pair)
models<- c("acg_naa", "acg_cr", "acg_cho", "acg_ins", "acg_glx", "lag_naa", "lag_cr", "lag_cho", "lag_ins", "lag_glx")

#list p-values for each model, matching order of models
p<-c(.0217, .307, .0258, .0494, .23, .00093, .115, .00229, .636, 0)

fdr_table<-data.frame(models, p)

fdr_table$p_fdr<-p.adjust(p, method = "fdr", n= length(p))
fdr_table
##     models       p       p_fdr
## 1  acg_naa 0.02170 0.051600000
## 2   acg_cr 0.30700 0.341111111
## 3  acg_cho 0.02580 0.051600000
## 4  acg_ins 0.04940 0.082333333
## 5  acg_glx 0.23000 0.287500000
## 6  lag_naa 0.00093 0.004650000
## 7   lag_cr 0.11500 0.164285714
## 8  lag_cho 0.00229 0.007633333
## 9  lag_ins 0.63600 0.636000000
## 10 lag_glx 0.00000 0.000000000
#run FDR without LAG-Glx in case it's invalid:
p2<-c(.0217, .307, .0258, .0494, .23, .00093, .115, .00229, .636)
p.adjust(p2, method = "fdr", n= length(p2))
## [1] 0.0580500 0.3453750 0.0580500 0.0889200 0.2957143 0.0083700 0.1725000
## [8] 0.0103050 0.6360000

##Additional models Next steps following OHBM abstract submission (Dec. 2021): test random slopes (LAG only)

#read in data
lag_dat<-read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320.csv")

#basic age model NAA
lag_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##      896      911     -444      888      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4313 -0.5792 -0.0185  0.5440  3.6263 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1641   0.4051  
##  Residual             0.8272   0.9095  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  13.82375    0.18588 309.72511  74.369  < 2e-16 ***
## mri_age_y     0.09699    0.02902 316.51634   3.343  0.00093 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.929
#test random slope NAA
lag_naa_age_randslope <- lmer(NAA_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id), data=lag_dat, REML=FALSE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00298042 (tol = 0.002, component 1)
summary(lag_naa_age_randslope)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    899.3    921.9   -443.7    887.3      312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4425 -0.5829 -0.0345  0.5457  3.6775 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subj_id  (Intercept) 0.438346 0.66208       
##           mri_age_y   0.005125 0.07159  -0.81
##  Residual             0.808385 0.89910       
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 13.81838    0.19532 84.82239  70.749  < 2e-16 ***
## mri_age_y    0.09843    0.03011 76.83548   3.269  0.00162 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.936
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00298042 (tol = 0.002, component 1)
#model failed to converge 

#basic age model Cre
lag_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    848.2    863.3   -420.1    840.2      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9764 -0.5307 -0.1261  0.5525  4.1548 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08477  0.2912  
##  Residual             0.73683  0.8584  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.49912    0.16805 308.49608  62.476   <2e-16 ***
## mri_age_y    -0.04183    0.02646 319.87112  -1.581    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.937
#test random slope Cre
lag_cre_age_randslope<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id), data=lag_dat, REML=FALSE)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00521307 (tol = 0.002, component 1)
summary(lag_cre_age_randslope)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    848.6    871.2   -418.3    836.6      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1240 -0.5211 -0.1430  0.5938  4.3226 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subj_id  (Intercept) 0.594776 0.77122       
##           mri_age_y   0.008986 0.09479  -0.96
##  Residual             0.705980 0.84023       
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 10.50247    0.18591 91.11070  56.492   <2e-16 ***
## mri_age_y   -0.04111    0.02825 82.73282  -1.455    0.149    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.951
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00521307 (tol = 0.002, component 1)
#model failed to converge

#basic age model Cho
lag_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     51.3     66.4    -21.7     43.3      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6539 -0.5810 -0.0562  0.4659  4.1202 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.0149   0.1221  
##  Residual             0.0561   0.2368  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.346670   0.049038 310.144986  47.854  < 2e-16 ***
## mri_age_y    -0.023315   0.007581 315.773363  -3.075  0.00229 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.921
#test random slope Cho
lag_cho_age_randslope<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_randslope)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     55.0     77.6    -21.5     43.0      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6693 -0.5804 -0.0496  0.4832  4.0955 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  subj_id  (Intercept) 0.0211331 0.14537       
##           mri_age_y   0.0002882 0.01698  -0.56
##  Residual             0.0550144 0.23455       
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  2.350155   0.050009 57.473557  46.994  < 2e-16 ***
## mri_age_y   -0.024037   0.007888 61.859198  -3.047  0.00339 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.924
#compare models
anova(lag_cho_age, lag_cho_age_randslope) #fit did not improve with inclusion of random slopes
## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_age_randslope: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##                       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_cho_age              4 51.345 66.419 -21.673   43.345                     
## lag_cho_age_randslope    6 55.014 77.624 -21.507   43.014 0.3314  2     0.8473
#basic age model Ins
lag_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1013.3   1028.3   -502.6   1005.3      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3984 -0.5594 -0.0111  0.5911  4.6516 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1768   0.4205  
##  Residual             1.2227   1.1057  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.28473    0.21938 307.28604  28.648   <2e-16 ***
## mri_age_y     0.01632    0.03444 318.45543   0.474    0.636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.934
#test random slope Ins
lag_ins_age_randslope<- lmer(Ins_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_ins_age_randslope)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1016.6   1039.1   -502.3   1004.6      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2826 -0.5584 -0.0263  0.5849  4.6042 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subj_id  (Intercept) 0.414651 0.64393       
##           mri_age_y   0.001291 0.03593  -1.00
##  Residual             1.211663 1.10076       
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.29900    0.22582 145.84527  27.894   <2e-16 ***
## mri_age_y     0.01355    0.03436 307.59889   0.394    0.694    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.937
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#boundary fit is singular

anova(lag_ins_age, lag_ins_age_randslope) #fit did not improve with inclusion of random slopes
## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_age_randslope: Ins_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##                       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_ins_age              4 1013.3 1028.3 -502.64   1005.3                     
## lag_ins_age_randslope    6 1016.6 1039.2 -502.28   1004.6 0.7233  2     0.6965
lag_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1396.1   1411.1   -694.0   1388.1      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6700 -0.6848  0.0127  0.6729  3.5338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             4.542    2.131   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  23.23676    0.39326 319.00000  59.088  < 2e-16 ***
## mri_age_y    -0.25961    0.06309 319.00000  -4.115 4.93e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#boundary fit is singular

lag_glx_age_randslope<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age_randslope)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1399.5   1422.1   -693.7   1387.5      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6292 -0.6832  0.0209  0.6274  3.5685 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subj_id  (Intercept) 1.32593  1.1515        
##           mri_age_y   0.03495  0.1869   -1.00
##  Residual             4.41427  2.1010        
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 23.19469    0.41423 77.21723  55.994  < 2e-16 ***
## mri_age_y   -0.24923    0.06668 67.41877  -3.738 0.000385 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.958
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#boundary fit is singular
anova(lag_glx_age, lag_glx_age_randslope) #fit did not improve with inclusion of random slopes
## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_age_randslope: Glu_Gln_conc_molal ~ mri_age_y + (1 + mri_age_y | subj_id)
##                       npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_glx_age              4 1396.1 1411.1 -694.04   1388.1                     
## lag_glx_age_randslope    6 1399.5 1422.1 -693.74   1387.5 0.5942  2      0.743

##Tissue fraction analysis Calculate tissue fraction GM/(GM+WM) for each subject/voxel Add tissue fraction as fixed effect in lmer models ACG

library(ggplot2)
library(EnvStats)
library(lme4)
library(lmerTest)
#ACG
#read in data
acg_dat <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113.csv")

#create tissue fraction variable
acg_dat$TissueFraction_GMWM <- acg_dat$fGM/(acg_dat$fGM + acg_dat$fWM)

#plot tissue fraction to check for issues
plot_ACG_tissuefraction <- ggplot(acg_dat, aes(TissueFraction_GMWM))
plot_ACG_tissuefraction+ geom_density() #big skew

plot_ACG_tissuefraction + geom_boxplot() #4 potential outliers

rosnerTest(acg_dat$TissueFraction_GMWM, k = 4)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4 
## 7.315531 5.204705 3.854675 3.544898 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 4 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 
## 3.425263 3.422302 3.419309 3.416284 
## 
## $n.outliers
## [1] 4
## 
## $alternative
## [1] "Up to 4 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 0.9390198 0.9661804 0.9655023 0.9366803 0.8741955 0.9311618 0.9320337
##   [8] 0.9397518 0.9548737 0.9669866 0.9444564 0.9138117 0.9490174 0.9573432
##  [15] 0.9497333 0.9748258 0.9941847 0.9486011 0.9602519 0.9612505 0.9732704
##  [22] 0.9790679 0.9479635 0.8046696 0.9719504 0.9435072 0.9860912 0.9252084
##  [29] 0.8946503 0.9275092 0.9413340 0.9619213 0.9344712 0.9996141 0.9611735
##  [36] 0.9812598 0.6594113 0.9283294 0.9853063 0.9976770 0.9481254 0.9628117
##  [43] 0.9347577 0.9444220 0.9509674 0.9248160 0.9468277 0.9750073 0.9871912
##  [50] 0.9848082 0.9777584 0.9709937 0.9427788 0.9516720 0.9772194 0.8991193
##  [57] 0.9724672 0.9533810 0.9848684 0.9562166 0.9898357 0.9339068 0.9541988
##  [64] 0.9779450 0.9755787 0.9424898 0.9860622 0.9440049 0.9238023 0.9511135
##  [71] 0.9519491 0.9895463 0.9500137 0.9847202 0.9509092 0.9630685 0.9736999
##  [78] 0.9281849 0.9266655 0.9872865 0.9798646 0.9595847 0.9541502 0.9831991
##  [85] 0.9377129 0.9309479 0.9847561 0.9729524 0.9754696 0.9392829 0.9671345
##  [92] 0.9700358 0.8595117 0.9559970 0.9768332 0.9787093 0.9556623 0.9457820
##  [99] 0.9936215 0.9881743 0.9804760 0.9764555 0.9591807 0.9290081 0.9073076
## [106] 0.9651388 0.9776632 0.9805040 0.9636278 0.9544220 0.9467918 0.9651846
## [113] 0.9882172
## 
## $data.name
## [1] "acg_dat$TissueFraction_GMWM"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i    Mean.i       SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 0.9530867 0.04014410 0.6594113      37 7.315531   3.425263    TRUE
## 2 1 0.9557088 0.02901974 0.8046696      24 5.204705   3.422302    TRUE
## 3 2 0.9570695 0.02530896 0.8595117      93 3.854675   3.419309    TRUE
## 4 3 0.9579564 0.02362857 0.8741955       5 3.544898   3.416284    TRUE
## 
## attr(,"class")
## [1] "gofOutlier"
#4 outliers: Obs. Num 37 (PS14_090), 24 (PS14_050), 93 (PS18_1601), 5 (PS14_009) - double check cogregistration
#PS14_090 is shifted to the left and includes substantial WM - EXCLUDE
#PS14_050 is shifted to the left and T1 quality is poor - EXCLUDE
#PS18_1601 head tilted so voxel captures a bit of WM in the left hemi., a little anterior; doesn't seem terrible
#PS14_009 looks like chin was tipped down toward chest a bit, corner of voxel captures some CC white matter but seems ok

#save database with tissue fraction
write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113_tissuefraction.csv", row.names = FALSE)

#exclude scans PS14_090 and PS14_050 due to outlier status coregistration issues
acg_dat_tf<-read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113_tissuefraction.csv")
acg_dat_tf["37", "TissueFraction_GMWM"] <- NA
acg_dat_tf["24", "TissueFraction_GMWM"] <- NA

acg_dat_tf <- subset(acg_dat_tf, study_code != c("PS14_090", "PS14_050"))
## Warning in study_code != c("PS14_090", "PS14_050"): longer object length is not
## a multiple of shorter object length
#check correlation between tissue fraction and age
acg_TisFrac_age <- lmer(TissueFraction_GMWM ~ mri_age_y + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_TisFrac_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: TissueFraction_GMWM ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##   -499.0   -488.1    253.5   -507.0      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.98848 -0.32600  0.01794  0.42686  1.85387 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  subj_id  (Intercept) 0.0004829 0.02197 
##  Residual             0.0001820 0.01349 
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 9.552e-01  8.892e-03 1.002e+02 107.423   <2e-16 ***
## mri_age_y   3.986e-04  2.077e-03 9.572e+01   0.192    0.848    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.957
#not sig.

#Add tissue fraction as fixed effect in ACG models
#NAA
#basic age model
acg_naa_age <- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_naa_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    364.2    375.0   -178.1    356.2      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78812 -0.44480  0.04041  0.41080  1.88431 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.8872   0.9419  
##  Residual             0.6396   0.7998  
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  14.0978     0.4398 108.8460  32.052   <2e-16 ***
## mri_age_y     0.2044     0.1031 106.6834   1.982   0.0501 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.960
#add tissue fraction
acg_naa_age_tf <- lmer(NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_naa_age_tf) #no sig. effect of tissue fraction
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    366.1    379.7   -178.1    356.1      106 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78485 -0.44445  0.02803  0.41427  1.84724 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.9017   0.9496  
##  Residual             0.6272   0.7919  
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)          13.1663     4.4412 110.4728   2.965  0.00371 **
## mri_age_y             0.2051     0.1030 106.1791   1.990  0.04915 * 
## TissueFraction_GMWM   0.9703     4.6197 109.9872   0.210  0.83402   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.093       
## TssFrc_GMWM -0.995 -0.002
#compare models
anova(acg_naa_age, acg_naa_age_tf) #model fit not sig. different
## Data: acg_dat_tf
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_age_tf: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_naa_age       4 364.15 374.99 -178.08   356.15                     
## acg_naa_age_tf    5 366.11 379.66 -178.06   356.11 0.0424  1     0.8369
#Cre
#basic age model
acg_cre_age <- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    324.5    335.4   -158.3    316.5      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.01085 -0.43475  0.01479  0.41135  1.81014 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.5452   0.7383  
##  Residual             0.5109   0.7147  
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  11.84959    0.36874 110.02636  32.136   <2e-16 ***
## mri_age_y     0.05989    0.08656 108.21331   0.692     0.49    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.961
#add tissue fraction
acg_cre_age_tf <- lmer(Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_cre_age_tf)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    325.8    339.4   -157.9    315.8      106 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.05786 -0.42131  0.01284  0.42273  1.79300 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.5271   0.7260  
##  Residual             0.5203   0.7213  
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          14.88776    3.69640 109.55054   4.028 0.000104 ***
## mri_age_y             0.05874    0.08634 108.44592   0.680 0.497715    
## TissueFraction_GMWM  -3.17067    3.84260 108.86252  -0.825 0.411099    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.100       
## TssFrc_GMWM -0.995  0.004
#Cho
#basic age model
acg_cho_age <- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##     72.2     83.0    -32.1     64.2      106 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2025 -0.6306 -0.1413  0.6299  2.4918 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.02935  0.1713  
##  Residual             0.07690  0.2773  
## Number of obs: 110, groups:  subj_id, 99
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.76889    0.11894 109.88524  23.281   <2e-16 ***
## mri_age_y    -0.06836    0.02799 109.96492  -2.442   0.0162 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.963
#add tissue fraction
acg_cho_age_tf <- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_cho_age_tf)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 |  
##     subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##     72.7     86.2    -31.4     62.7      105 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2564 -0.5970 -0.1487  0.6628  2.5884 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.02784  0.1668  
##  Residual             0.07690  0.2773  
## Number of obs: 110, groups:  subj_id, 99
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           4.19678    1.17742 109.25030   3.564 0.000542 ***
## mri_age_y            -0.06876    0.02781 109.97149  -2.473 0.014935 *  
## TissueFraction_GMWM  -1.49055    1.22297 109.07759  -1.219 0.225551    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.109       
## TssFrc_GMWM -0.995  0.013
#compare models
anova(acg_cho_age, acg_cho_age_tf) #fit not significantly different
## Data: acg_dat_tf
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_age_tf: Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## acg_cho_age       4 72.191 82.993 -32.096   64.191                    
## acg_cho_age_tf    5 72.718 86.221 -31.359   62.718 1.473  1     0.2249
#Ins
#basic age model
acg_ins_age <- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    407.8    418.6   -199.9    399.8      106 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5246 -0.4721  0.0711  0.5130  1.5756 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.9237   0.9611  
##  Residual             1.3525   1.1630  
## Number of obs: 110, groups:  subj_id, 99
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   8.5924     0.5484 109.9754  15.667   <2e-16 ***
## mri_age_y    -0.2258     0.1294 109.4869  -1.746   0.0837 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.962
#add tissue fraction
acg_ins_age_tf <- lmer(Ins_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_ins_age_tf)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    409.4    422.9   -199.7    399.4      105 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.52269 -0.47942  0.08856  0.57698  1.49140 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.9414   0.9702  
##  Residual             1.3276   1.1522  
## Number of obs: 110, groups:  subj_id, 99
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)           4.8246     5.4937 108.8944   0.878   0.3818  
## mri_age_y            -0.2249     0.1291 109.3760  -1.742   0.0844 .
## TissueFraction_GMWM   3.9350     5.7034 108.5006   0.690   0.4917  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.116       
## TssFrc_GMWM -0.995  0.021
#Glx
#basic age model
acg_glx_age <- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    530.3    541.1   -261.1    522.3      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.48912 -0.37848 -0.00398  0.44809  1.68844 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.601    2.145   
##  Residual             2.344    1.531   
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  25.9373     0.9245 105.2183  28.055   <2e-16 ***
## mri_age_y     0.2288     0.2164 101.8687   1.057    0.293    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.959
#add tissue fraction
acg_glx_age_tf <- lmer(Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=acg_dat_tf, REML=FALSE)
summary(acg_glx_age_tf)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: acg_dat_tf
## 
##      AIC      BIC   logLik deviance df.resid 
##    530.5    544.0   -260.2    520.5      106 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.43181 -0.34955 -0.03061  0.41546  1.49601 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.837    2.199   
##  Residual             2.068    1.438   
## Number of obs: 111, groups:  subj_id, 100
## 
## Fixed effects:
##                     Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)          13.2424     9.3236 110.9911   1.420    0.158
## mri_age_y             0.2289     0.2135  98.4861   1.072    0.286
## TissueFraction_GMWM  13.2676     9.7105 110.9341   1.366    0.175
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.080       
## TssFrc_GMWM -0.995 -0.014

LAG Tissue Fraction models

#LAG
#read in data
lag_dat<-read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320.csv")

#create tissue fraction variable
lag_dat$TissueFraction_GMWM <- lag_dat$fGM/(lag_dat$fGM + lag_dat$fWM)

#plot tissue fraction to check for issues
plot_LAG_tissuefraction <- ggplot(lag_dat, aes(TissueFraction_GMWM))
plot_LAG_tissuefraction+ geom_density() #normal

plot_LAG_tissuefraction + geom_boxplot() #6 potential outliers

rosnerTest(lag_dat$TissueFraction_GMWM, k = 6) #no significant outliers
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6 
## 3.280842 3.075086 2.985173 2.935103 2.963794 2.932421 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 6 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 6 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 0.7657305 0.6515751 0.5744539 0.6894155 0.8381194 0.5951315 0.7008113
##   [8] 0.6813746 0.5352570 0.6558286 0.7097519 0.7819635 0.7021001 0.7039674
##  [15] 0.4764324 0.6028844 0.6327450 0.6087942 0.6968812 0.8026664 0.7009372
##  [22] 0.7208761 0.6605230 0.7324027 0.5453879 0.7620629 0.6701616 0.6375845
##  [29] 0.5632865 0.6708753 0.6927008 0.7798970 0.7618192 0.6477681 0.6613518
##  [36] 0.5340622 0.5942859 0.6943742 0.7009750 0.5879671 0.4411657 0.6156132
##  [43] 0.6142354 0.7259517 0.5979064 0.5685179 0.5271897 0.5355836 0.7237169
##  [50] 0.7571831 0.5285119 0.6751229 0.5119628 0.7521085 0.6615893 0.5425676
##  [57] 0.6627217 0.6549878 0.7279301 0.7256598 0.7477565 0.6395054 0.7213695
##  [64] 0.7384663 0.7633997 0.7529087 0.8002460 0.8145222 0.7366533 0.7977780
##  [71] 0.6151766 0.7154834 0.7754343 0.8142228 0.7207045 0.6896770 0.5855033
##  [78] 0.4519376 0.3757102 0.6561017 0.7878613 0.7377081 0.8047878 0.7816798
##  [85] 0.8279033 0.6753913 0.7666730 0.5795116 0.7313181 0.6599876 0.7201454
##  [92] 0.6696166 0.5834821 0.6485976 0.5740400 0.7069867 0.7142496 0.6394643
##  [99] 0.7231628 0.6282001 0.7420446 0.5742007 0.7048192 0.8067634 0.8402937
## [106] 0.6709200 0.7285484 0.6456249 0.7123905 0.7182381 0.6336725 0.7074584
## [113] 0.5114034 0.7548622 0.7512744 0.6852121 0.6741392 0.4619200 0.5944252
## [120] 0.6890812 0.5371624 0.6516809 0.6145975 0.6753211 0.5748713 0.7172832
## [127] 0.7275696 0.4961396 0.6162007 0.4003669 0.6812224 0.6283159 0.6387185
## [134] 0.7522739 0.4447576 0.5459764 0.5522673 0.5951632 0.5682089 0.5371529
## [141] 0.6615547 0.6791234 0.7036370 0.6886679 0.6631862 0.6341630 0.6738518
## [148] 0.7549798 0.6955773 0.6542445 0.6925320 0.6505739 0.6425831 0.6677008
## [155] 0.7839529 0.6553603 0.5432005 0.6122120 0.6548850 0.6367154 0.7679106
## [162] 0.7485923 0.7184663 0.6312047 0.7349290 0.7043475 0.4988189 0.8401854
## [169] 0.6577281 0.8438765 0.8390803 0.6282136 0.6405363 0.7191654 0.6334385
## [176] 0.7081752 0.6659769 0.6107064 0.6347357 0.7290806 0.7344315 0.6660147
## [183] 0.6251258 0.5997631 0.5785047 0.7984798 0.5070360 0.6098369 0.6299690
## [190] 0.6516707 0.5494253 0.6308773 0.7664214 0.5915382 0.5878598 0.8187175
## [197] 0.6741752 0.7564550 0.6443908 0.7395504 0.5204329 0.6170451 0.7002564
## [204] 0.6246470 0.6681476 0.5771348 0.7295070 0.7423764 0.6384295 0.5430082
## [211] 0.6707773 0.6881614 0.7608897 0.5543782 0.6804216 0.5896538 0.6427618
## [218] 0.6536086 0.7261262 0.6317773 0.6652328 0.6605223 0.5295042 0.6174694
## [225] 0.6425888 0.5149187 0.7329097 0.5653457 0.5104815 0.6635529 0.7475020
## [232] 0.7604907 0.6907149 0.5576525 0.6882708 0.7095940 0.5975811 0.6217704
## [239] 0.6017863 0.6938276 0.5893048 0.4639501 0.7554426 0.5394351 0.5730658
## [246] 0.5198421 0.6217217 0.5560520 0.7715035 0.7450276 0.7578195 0.7198068
## [253] 0.6385403 0.4769133 0.6216745 0.7215428 0.7135267 0.5015701 0.5463705
## [260] 0.3420581 0.4840778 0.4052625 0.3302661 0.3499400 0.6039142 0.4129804
## [267] 0.5315517 0.5809198 0.4792521 0.5046752 0.5684662 0.5218814 0.3154764
## [274] 0.3403429 0.5024469 0.3751392 0.6127580 0.7315085 0.4137321 0.4513800
## [281] 0.5718875 0.7465811 0.5551458 0.4922670 0.5978026 0.6036407 0.7885750
## [288] 0.7143374 0.5734673 0.7879821 0.5626092 0.6622175 0.4524319 0.4110109
## [295] 0.6597638 0.5684290 0.4234816 0.4456416 0.5029994 0.5793653 0.6290320
## [302] 0.4115474 0.6699743 0.7957488 0.7141763 0.7477478 0.5862060 0.7635972
## [309] 0.7800932 0.5786560 0.5839172 0.6023700 0.5311156 0.5390859 0.6565085
## [316] 0.8540765 0.7585903 0.5580707 0.6451797 0.4621925 0.9951620
## 
## $data.name
## [1] "lag_dat$TissueFraction_GMWM"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i    Mean.i      SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 0.6422445 0.1075692 0.9951620     321 3.280842   3.742579   FALSE
## 2 1 0.6411416 0.1059044 0.3154764     273 3.075086   3.741706   FALSE
## 3 2 0.6421625 0.1044818 0.3302661     263 2.985173   3.740829   FALSE
## 4 3 0.6431433 0.1031652 0.3403429     274 2.935103   3.739949   FALSE
## 5 4 0.6440985 0.1019101 0.3420581     260 2.963794   3.739067   FALSE
## 6 5 0.6450543 0.1006385 0.3499400     264 2.932421   3.738181   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320_tissuefraction.csv")

#check correlation between tissue fraction and age
lag_TisFrac_age <- lmer(TissueFraction_GMWM ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_TisFrac_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: TissueFraction_GMWM ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   -562.8   -547.7    285.4   -570.8      317 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7967 -0.6352  0.0458  0.6787  3.5935 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  subj_id  (Intercept) 0.0006557 0.02561 
##  Residual             0.0093041 0.09646 
## Number of obs: 321, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   0.770798   0.018475 306.068797  41.722  < 2e-16 ***
## mri_age_y    -0.021705   0.002927 320.934582  -7.417 1.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.943
#tissue fraction is significantly negatively associated with age

#in a simple correlation test
cor.test(lag_dat$mri_age_y, lag_dat$TissueFraction_GMWM) # also significant: r = -.37
## 
##  Pearson's product-moment correlation
## 
## data:  lag_dat$mri_age_y and lag_dat$TissueFraction_GMWM
## t = -7.1217, df = 319, p-value = 7.113e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4611519 -0.2719374
## sample estimates:
##        cor 
## -0.3703805
#Add tissue fraction as a fixed effect in LAG models
#basic age model NAA
lag_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)

#add tissue fraction
lag_naa_age_TisFrac <- lmer(NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age_TisFrac)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    858.0    876.8   -424.0    848.0      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9925 -0.5718 -0.0199  0.5594  3.9440 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1630   0.4037  
##  Residual             0.7184   0.8476  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          11.22872    0.43459 317.97211  25.838  < 2e-16 ***
## mri_age_y             0.17233    0.02961 317.97750   5.819 1.44e-08 ***
## TissueFraction_GMWM   3.35162    0.51255 313.38178   6.539 2.51e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.703       
## TssFrc_GMWM -0.915  0.394
#compare model fit
anova(lag_naa_age, lag_naa_age_TisFrac) #fit significantly improved
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age_TisFrac: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_naa_age            4 895.96 911.01 -443.98   887.96                     
## lag_naa_age_TisFrac    5 857.96 876.77 -423.98   847.96 39.996  1  2.545e-10
##                        
## lag_naa_age            
## lag_naa_age_TisFrac ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add tissue fraction interaction with age
lag_naa_age_TisFracXage <- lmer(NAA_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age_TisFracXage)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    857.3    879.9   -422.6    845.3      312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1567 -0.5361 -0.0447  0.5566  4.0061 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1599   0.3999  
##  Residual             0.7134   0.8447  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##                               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                     9.7907     0.9790 313.2768  10.001  < 2e-16 ***
## mri_age_y                       0.3877     0.1348 307.1007   2.877 0.004302 ** 
## TissueFraction_GMWM             5.7025     1.5239 311.3849   3.742 0.000217 ***
## mri_age_y:TissueFraction_GMWM  -0.3606     0.2202 305.3472  -1.638 0.102532    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.943              
## TssFrc_GMWM -0.981  0.948       
## m__:TF_GMWM  0.897 -0.976 -0.942
#interaction term for agexTissueFraction not sig.

#compare model fit
anova(lag_naa_age_TisFrac, lag_naa_age_TisFracXage) #not sig., stick with age + TissueFraction model
## Data: lag_dat
## Models:
## lag_naa_age_TisFrac: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_naa_age_TisFracXage: NAA_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_naa_age_TisFrac        5 857.96 876.77 -423.98   847.96          
## lag_naa_age_TisFracXage    6 857.29 879.86 -422.65   845.29 2.6694  1
##                         Pr(>Chisq)
## lag_naa_age_TisFrac               
## lag_naa_age_TisFracXage     0.1023
#add sex to tissue fraction model
lag_naa_age_TisFrac_sex <- lmer(NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age_TisFrac_sex) #sex not significant
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    859.9    882.4   -423.9    847.9      312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9834 -0.5775 -0.0126  0.5684  3.9313 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1616   0.4020  
##  Residual             0.7190   0.8479  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          11.24098    0.43642 317.96672  25.757  < 2e-16 ***
## mri_age_y             0.17281    0.02965 317.99247   5.828 1.37e-08 ***
## TissueFraction_GMWM   3.35735    0.51294 313.49233   6.545 2.42e-10 ***
## female               -0.03974    0.13192  97.72196  -0.301    0.764    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.695              
## TssFrc_GMWM -0.907  0.395       
## female      -0.093 -0.052 -0.039
#compare model fit
anova(lag_naa_age_TisFrac, lag_naa_age_TisFrac_sex) #not sig., stick with age + TissueFraction model
## Data: lag_dat
## Models:
## lag_naa_age_TisFrac: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_naa_age_TisFrac_sex: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_naa_age_TisFrac        5 857.96 876.77 -423.98   847.96          
## lag_naa_age_TisFrac_sex    6 859.87 882.44 -423.94   847.87 0.0901  1
##                         Pr(>Chisq)
## lag_naa_age_TisFrac               
## lag_naa_age_TisFrac_sex     0.7641
#Cre
#basic age model
lag_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    848.2    863.3   -420.1    840.2      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9764 -0.5307 -0.1261  0.5525  4.1548 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08477  0.2912  
##  Residual             0.73683  0.8584  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.49912    0.16805 308.49608  62.476   <2e-16 ***
## mri_age_y    -0.04183    0.02646 319.87112  -1.581    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.937
#add tissue fraction
lag_cre_age_TisFrac<- lmer(Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age_TisFrac)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    760.2    779.0   -375.1    750.2      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7493 -0.5056  0.0270  0.4419  4.8839 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08547  0.2923  
##  Residual             0.54120  0.7357  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           7.04677    0.36923 319.70413  19.085   <2e-16 ***
## mri_age_y             0.05661    0.02497 319.62390   2.267   0.0241 *  
## TissueFraction_GMWM   4.47115    0.43745 317.80861  10.221   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.701       
## TssFrc_GMWM -0.918  0.393
#compare model fit
anova(lag_cre_age, lag_cre_age_TisFrac) #sig.
## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_age_TisFrac: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_cre_age            4 848.23 863.31 -420.12   840.23                     
## lag_cre_age_TisFrac    5 760.20 779.04 -375.10   750.20 90.029  1  < 2.2e-16
##                        
## lag_cre_age            
## lag_cre_age_TisFrac ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add interaction
lag_cre_age_TisFrac_int<- lmer(Cre_PCr_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age_TisFrac_int)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    761.9    784.5   -375.0    749.9      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7360 -0.5018  0.0276  0.4507  4.9009 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08463  0.2909  
##  Residual             0.54125  0.7357  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                                Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                     6.65386    0.83511 317.57482   7.968 2.93e-14
## mri_age_y                       0.11569    0.11532 312.49738   1.003 0.316524
## TissueFraction_GMWM             5.11282    1.29936 316.10788   3.935 0.000102
## mri_age_y:TissueFraction_GMWM  -0.09884    0.18836 310.83144  -0.525 0.600131
##                                  
## (Intercept)                   ***
## mri_age_y                        
## TissueFraction_GMWM           ***
## mri_age_y:TissueFraction_GMWM    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.943              
## TssFrc_GMWM -0.981  0.948       
## m__:TF_GMWM  0.897 -0.976 -0.942
#compare model fit
anova(lag_cre_age_TisFrac, lag_cre_age_TisFrac_int) #not sig.
## Data: lag_dat
## Models:
## lag_cre_age_TisFrac: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_cre_age_TisFrac_int: Cre_PCr_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_cre_age_TisFrac        5 760.20 779.04 -375.10   750.20          
## lag_cre_age_TisFrac_int    6 761.93 784.54 -374.96   749.93 0.2748  1
##                         Pr(>Chisq)
## lag_cre_age_TisFrac               
## lag_cre_age_TisFrac_int     0.6002
#add sex to tissue fraction model

lag_cre_age_TisFrac_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age_TisFrac_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    762.0    784.6   -375.0    750.0      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7681 -0.5123  0.0220  0.4570  4.9061 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.0852   0.2919  
##  Residual             0.5410   0.7355  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           7.03303    0.37047 319.29640  18.984   <2e-16 ***
## mri_age_y             0.05601    0.02501 319.56737   2.240   0.0258 *  
## TissueFraction_GMWM   4.46392    0.43763 317.89203  10.200   <2e-16 ***
## female                0.04644    0.10629  90.75521   0.437   0.6632    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.693              
## TssFrc_GMWM -0.911  0.394       
## female      -0.085 -0.056 -0.037
#compare model fit
anova(lag_cre_age_TisFrac, lag_cre_age_TisFrac_sex) #not sig.
## Data: lag_dat
## Models:
## lag_cre_age_TisFrac: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_cre_age_TisFrac_sex: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_cre_age_TisFrac        5 760.20 779.04 -375.10   750.20          
## lag_cre_age_TisFrac_sex    6 762.01 784.62 -375.01   750.01 0.1908  1
##                         Pr(>Chisq)
## lag_cre_age_TisFrac               
## lag_cre_age_TisFrac_sex     0.6622
#Cho
#basic age model
lag_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     51.3     66.4    -21.7     43.3      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6539 -0.5810 -0.0562  0.4659  4.1202 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.0149   0.1221  
##  Residual             0.0561   0.2368  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.346670   0.049038 310.144986  47.854  < 2e-16 ***
## mri_age_y    -0.023315   0.007581 315.773363  -3.075  0.00229 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.921
#add tissue fraction
lag_cho_age_TisFrac<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_TisFrac)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     25.3     44.2     -7.7     15.3      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5688 -0.6208 -0.0348  0.5176  4.6606 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01202  0.1096  
##  Residual             0.05232  0.2287  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           2.929074   0.117156 319.933789  25.002  < 2e-16 ***
## mri_age_y            -0.040403   0.007925 319.967590  -5.098 5.89e-07 ***
## TissueFraction_GMWM  -0.749693   0.138214 313.218002  -5.424 1.17e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.704       
## TssFrc_GMWM -0.917  0.399
#compare model fit
anova(lag_cho_age, lag_cho_age_TisFrac) #sig
## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_age_TisFrac: Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                     npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)
## lag_cho_age            4 51.345 66.419 -21.6727   43.345                     
## lag_cho_age_TisFrac    5 25.333 44.174  -7.6663   15.333 28.013  1  1.205e-07
##                        
## lag_cho_age            
## lag_cho_age_TisFrac ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add interaction
lag_cho_age_TisFrac_int<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_TisFrac_int) #interaction is sig. (p=.046)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     23.3     46.0     -5.7     11.3      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5480 -0.6259 -0.0670  0.5067  4.4339 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01168  0.1081  
##  Residual             0.05179  0.2276  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                                Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                     2.45851    0.26231 312.15748   9.373   <2e-16
## mri_age_y                       0.03027    0.03614 303.65993   0.838    0.403
## TissueFraction_GMWM             0.01894    0.40778 309.31482   0.046    0.963
## mri_age_y:TissueFraction_GMWM  -0.11820    0.05899 301.29296  -2.004    0.046
##                                  
## (Intercept)                   ***
## mri_age_y                        
## TissueFraction_GMWM              
## mri_age_y:TissueFraction_GMWM *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.943              
## TssFrc_GMWM -0.981  0.948       
## m__:TF_GMWM  0.896 -0.976 -0.942
#compare model fit
anova(lag_cho_age_TisFrac, lag_cho_age_TisFrac_int) #sig. (p<.046)
## Data: lag_dat
## Models:
## lag_cho_age_TisFrac: Cho_GPC_PCh_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_cho_age_TisFrac_int: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_cho_age_TisFrac        5 25.333 44.174 -7.6663   15.333          
## lag_cho_age_TisFrac_int    6 23.345 45.955 -5.6724   11.345 3.9878  1
##                         Pr(>Chisq)  
## lag_cho_age_TisFrac                 
## lag_cho_age_TisFrac_int    0.04583 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add sex to model
lag_cho_age_TisFrac_int_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_TisFrac_int_sex) #not sig.
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     24.2     50.6     -5.1     10.2      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5865 -0.6287 -0.0611  0.5043  4.4825 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01157  0.1076  
##  Residual             0.05164  0.2273  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                                 Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                     2.485820   0.263137 312.684444   9.447   <2e-16
## mri_age_y                       0.028390   0.036124 303.303507   0.786   0.4325
## TissueFraction_GMWM            -0.001475   0.407568 309.052784  -0.004   0.9971
## female                         -0.037481   0.035376  77.480663  -1.060   0.2927
## mri_age_y:TissueFraction_GMWM  -0.114251   0.059018 300.819739  -1.936   0.0538
##                                  
## (Intercept)                   ***
## mri_age_y                        
## TissueFraction_GMWM              
## female                           
## mri_age_y:TissueFraction_GMWM .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW female
## mri_age_y   -0.942                     
## TssFrc_GMWM -0.980  0.948              
## female      -0.098  0.050  0.047       
## m__:TF_GMWM  0.896 -0.976 -0.942 -0.064
#compare model fit
anova(lag_cho_age_TisFrac_int, lag_cho_age_TisFrac_int_sex) #not sig.
## Data: lag_dat
## Models:
## lag_cho_age_TisFrac_int: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
## lag_cho_age_TisFrac_int_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + female + (1 | subj_id)
##                             npar    AIC    BIC  logLik deviance  Chisq Df
## lag_cho_age_TisFrac_int        6 23.345 45.955 -5.6724   11.345          
## lag_cho_age_TisFrac_int_sex    7 24.224 50.603 -5.1122   10.225 1.1204  1
##                             Pr(>Chisq)
## lag_cho_age_TisFrac_int               
## lag_cho_age_TisFrac_int_sex     0.2898
#Ins
#basic age model
lag_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1013.3   1028.3   -502.6   1005.3      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3984 -0.5594 -0.0111  0.5911  4.6516 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1768   0.4205  
##  Residual             1.2227   1.1057  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.28473    0.21938 307.28604  28.648   <2e-16 ***
## mri_age_y     0.01632    0.03444 318.45543   0.474    0.636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.934
#add tissue fraction
lag_ins_age_TisFrac<- lmer(Ins_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age_TisFrac) #tissue fraction sig.; age not sig.
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1005.5   1024.3   -497.7    995.5      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3841 -0.5608 -0.0382  0.5912  4.7226 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1568   0.396   
##  Residual             1.1958   1.094   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           4.70590    0.54440 318.16597   8.644 2.69e-16 ***
## mri_age_y             0.06153    0.03685 318.01860   1.670  0.09596 .  
## TissueFraction_GMWM   2.04247    0.64587 317.48996   3.162  0.00172 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.700       
## TssFrc_GMWM -0.918  0.391
#compare model fit
anova(lag_ins_age, lag_ins_age_TisFrac) #sig.
## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_age_TisFrac: Ins_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                     npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)   
## lag_ins_age            4 1013.3 1028.3 -502.64  1005.28                       
## lag_ins_age_TisFrac    5 1005.5 1024.3 -497.73   995.47 9.812  1   0.001734 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add interaction

lag_ins_age_TisFrac_int<- lmer(Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age_TisFrac_int) #tissue fraction, age & interaction all sig.
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1001.6   1024.2   -494.8    989.6      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3911 -0.5979 -0.0053  0.5734  4.7179 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1767   0.4203  
##  Residual             1.1584   1.0763  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                     2.0286     1.2206 316.2405   1.662 0.097514 .  
## mri_age_y                       0.4654     0.1686 310.2011   2.761 0.006114 ** 
## TissueFraction_GMWM             6.4126     1.8995 314.4820   3.376 0.000828 ***
## mri_age_y:TissueFraction_GMWM  -0.6751     0.2755 308.0776  -2.450 0.014822 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.943              
## TssFrc_GMWM -0.981  0.948       
## m__:TF_GMWM  0.897 -0.976 -0.942
#compare model fit
anova(lag_ins_age_TisFrac, lag_ins_age_TisFrac_int) #sig.
## Data: lag_dat
## Models:
## lag_ins_age_TisFrac: Ins_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_ins_age_TisFrac_int: Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_ins_age_TisFrac        5 1005.5 1024.3 -497.73   995.47          
## lag_ins_age_TisFrac_int    6 1001.6 1024.2 -494.80   989.60 5.8652  1
##                         Pr(>Chisq)  
## lag_ins_age_TisFrac                 
## lag_ins_age_TisFrac_int    0.01544 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add sex
lag_ins_age_TisFrac_int_sex<- lmer(Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age_TisFrac_int_sex) #tissue fraction, age & interaction all sig.
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + female + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1003.6   1030.0   -494.8    989.6      312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3882 -0.5982 -0.0079  0.5748  4.7194 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1772   0.4209  
##  Residual             1.1581   1.0761  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                                 Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                     2.023474   1.226394 316.309945   1.650 0.099948
## mri_age_y                       0.465876   0.168831 309.440740   2.759 0.006135
## TissueFraction_GMWM             6.417150   1.901892 314.005179   3.374 0.000833
## female                          0.005839   0.155173  78.119816   0.038 0.970081
## mri_age_y:TissueFraction_GMWM  -0.675931   0.276096 307.149063  -2.448 0.014918
##                                  
## (Intercept)                   .  
## mri_age_y                     ** 
## TissueFraction_GMWM           ***
## female                           
## mri_age_y:TissueFraction_GMWM *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW female
## mri_age_y   -0.942                     
## TssFrc_GMWM -0.980  0.948              
## female      -0.097  0.054  0.051       
## m__:TF_GMWM  0.897 -0.976 -0.942 -0.068
#compare model fit
anova(lag_ins_age_TisFrac_int, lag_ins_age_TisFrac_int_sex) #not sig.
## Data: lag_dat
## Models:
## lag_ins_age_TisFrac_int: Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
## lag_ins_age_TisFrac_int_sex: Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + female + (1 | subj_id)
##                             npar    AIC    BIC logLik deviance  Chisq Df
## lag_ins_age_TisFrac_int        6 1001.6 1024.2 -494.8    989.6          
## lag_ins_age_TisFrac_int_sex    7 1003.6 1030.0 -494.8    989.6 0.0014  1
##                             Pr(>Chisq)
## lag_ins_age_TisFrac_int               
## lag_ins_age_TisFrac_int_sex     0.9704
#Glx
#basic age model
lag_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1396.1   1411.1   -694.0   1388.1      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6700 -0.6848  0.0127  0.6729  3.5338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             4.542    2.131   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  23.23676    0.39326 319.00000  59.088  < 2e-16 ***
## mri_age_y    -0.25961    0.06309 319.00000  -4.115 4.93e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#add tissue fraction
lag_glx_age_TisFrac<- lmer(Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age_TisFrac)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1238.1   1257.0   -614.1   1228.1      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5542 -0.6325 -0.0043  0.6145  4.7331 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             2.751    1.659   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          12.97325    0.77521 319.00000  16.735   <2e-16 ***
## mri_age_y             0.02096    0.05282 319.00000   0.397    0.692    
## TissueFraction_GMWM  13.37772    0.92834 319.00000  14.410   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.688       
## TssFrc_GMWM -0.919  0.369
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#compare models
anova(lag_glx_age, lag_glx_age_TisFrac) #sig
## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_age_TisFrac: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_glx_age            4 1396.1 1411.1 -694.04   1388.1                     
## lag_glx_age_TisFrac    5 1238.1 1257.0 -614.07   1228.1 159.94  1  < 2.2e-16
##                        
## lag_glx_age            
## lag_glx_age_TisFrac ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#add interaction
lag_glx_age_TisFrac_int<- lmer(Glu_Gln_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age_TisFrac_int)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1237.7   1260.3   -612.8   1225.7      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7571 -0.6187  0.0099  0.5924  4.7144 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.00     0.000   
##  Residual             2.73     1.652   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                    15.4568     1.7652 319.0000   8.756  < 2e-16 ***
## mri_age_y                      -0.3547     0.2458 319.0000  -1.443 0.149993    
## TissueFraction_GMWM             9.3213     2.7527 319.0000   3.386 0.000797 ***
## mri_age_y:TissueFraction_GMWM   0.6292     0.4022 319.0000   1.565 0.118677    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.943              
## TssFrc_GMWM -0.982  0.947       
## m__:TF_GMWM  0.899 -0.977 -0.942
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#compare models
anova(lag_glx_age_TisFrac, lag_glx_age_TisFrac_int) #not sig
## Data: lag_dat
## Models:
## lag_glx_age_TisFrac: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_glx_age_TisFrac_int: Glu_Gln_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_glx_age_TisFrac        5 1238.1 1257.0 -614.07   1228.1          
## lag_glx_age_TisFrac_int    6 1237.7 1260.3 -612.85   1225.7 2.4385  1
##                         Pr(>Chisq)
## lag_glx_age_TisFrac               
## lag_glx_age_TisFrac_int     0.1184
#add sex
lag_glx_age_TisFrac_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age_TisFrac_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1239.1   1261.7   -613.5   1227.1      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4999 -0.6414  0.0131  0.6445  4.7946 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             2.742    1.656   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          12.91764    0.77583 319.00000  16.650   <2e-16 ***
## mri_age_y             0.01866    0.05278 319.00000   0.353    0.724    
## TissueFraction_GMWM  13.34681    0.92730 319.00000  14.393   <2e-16 ***
## female                0.19087    0.18604 319.00000   1.026    0.306    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.683              
## TssFrc_GMWM -0.914  0.369       
## female      -0.070 -0.042 -0.032
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
#compare models
anova(lag_glx_age_TisFrac, lag_glx_age_TisFrac_sex) #not sig
## Data: lag_dat
## Models:
## lag_glx_age_TisFrac: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
## lag_glx_age_TisFrac_sex: Glu_Gln_conc_molal ~ mri_age_y + TissueFraction_GMWM + female + (1 | subj_id)
##                         npar    AIC    BIC  logLik deviance  Chisq Df
## lag_glx_age_TisFrac        5 1238.1 1257.0 -614.07   1228.1          
## lag_glx_age_TisFrac_sex    6 1239.1 1261.7 -613.54   1227.1 1.0509  1
##                         Pr(>Chisq)
## lag_glx_age_TisFrac               
## lag_glx_age_TisFrac_sex     0.3053
#check models for multicollinearity
library(performance)

check_collinearity(lag_naa_age_TisFrac) #Low
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                 Term  VIF Increased SE Tolerance
##            mri_age_y 1.18         1.09      0.84
##  TissueFraction_GMWM 1.18         1.09      0.84
check_collinearity(lag_cre_age_TisFrac) #Low
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                 Term  VIF Increased SE Tolerance
##            mri_age_y 1.18         1.09      0.85
##  TissueFraction_GMWM 1.18         1.09      0.85
check_collinearity(lag_cho_age_TisFrac_int) #interaction term causes inflated VIFs - check model without interaction
## Warning: Model has interaction terms. VIFs might be inflated. You may check
##   multicollinearity among predictors of a model without interaction terms.
## # Check for Multicollinearity
## 
## High Correlation
## 
##                           Term   VIF Increased SE Tolerance
##                      mri_age_y 25.02         5.00      0.04
##            TissueFraction_GMWM 10.47         3.24      0.10
##  mri_age_y:TissueFraction_GMWM 22.26         4.72      0.04
check_collinearity(lag_cho_age_TisFrac) #Low
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                 Term  VIF Increased SE Tolerance
##            mri_age_y 1.19         1.09      0.84
##  TissueFraction_GMWM 1.19         1.09      0.84
check_collinearity(lag_ins_age_TisFrac) #Low
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                 Term  VIF Increased SE Tolerance
##            mri_age_y 1.18         1.09      0.85
##  TissueFraction_GMWM 1.18         1.09      0.85
check_collinearity(lag_glx_age_TisFrac) #Low
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                 Term  VIF Increased SE Tolerance
##            mri_age_y 1.16         1.08      0.86
##  TissueFraction_GMWM 1.16         1.08      0.86

##Unpack interactions Use sim_slopes function from interactions package to conduct simple slopes analysis on continuous x continuous interactions

library(interactions)

#LAG Cho: age x tissue fraction interaction
lag_cho_age_TisFrac_int<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)

#run simple slopes analysis with Johnson-Neyman plot
sim_slopes(lag_cho_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM, jnplot = TRUE)
## JOHNSON-NEYMAN INTERVAL 
## 
## When TissueFraction_GMWM is OUTSIDE the interval [-18.37, 0.45], the slope
## of mri_age_y is p < .05.
## 
## Note: The range of observed values of TissueFraction_GMWM is [0.32, 1.00]

## SIMPLE SLOPES ANALYSIS 
## 
## Slope of mri_age_y when TissueFraction_GMWM = 0.5347131 (- 1 SD): 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.03   0.01    -3.78   0.00
## 
## Slope of mri_age_y when TissueFraction_GMWM = 0.6424099 (Mean): 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.05   0.01    -5.50   0.00
## 
## Slope of mri_age_y when TissueFraction_GMWM = 0.7501066 (+ 1 SD): 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.06   0.01    -4.89   0.00
#save Johnson-Neyman plot
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_cho_age_tf_int_JNplot.png")
## Saving 7 x 5 in image
#create interaction plot (make it pretty later)
#basic plot
Cho_age_tf_int_plot <- interact_plot(lag_cho_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM)
Cho_age_tf_int_plot

#linearity check
Cho_age_tf_int_plot_lin <- interact_plot(lag_cho_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM, modx.values="terciles", plot.points=TRUE, linearity.check=TRUE)
## Medians of each tercile of TissueFraction_GMWM are 0.539, 0.655, 0.747
Cho_age_tf_int_plot_lin
## `geom_smooth()` using formula 'y ~ x'

#LAG Ins: age x tissue fraction interaction
lag_ins_age_TisFrac_int<- lmer(Ins_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 | subj_id), data=lag_dat, REML=FALSE)

#run simple slopes analysis with Johnson-Neyman plot
sim_slopes(lag_ins_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM, jnplot = TRUE)
## JOHNSON-NEYMAN INTERVAL 
## 
## When TissueFraction_GMWM is OUTSIDE the interval [0.58, 1.13], the slope of
## mri_age_y is p < .05.
## 
## Note: The range of observed values of TissueFraction_GMWM is [0.32, 1.00]

## SIMPLE SLOPES ANALYSIS 
## 
## Slope of mri_age_y when TissueFraction_GMWM = 0.5344591 (- 1 SD): 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   0.10   0.04     2.59   0.01
## 
## Slope of mri_age_y when TissueFraction_GMWM = 0.6423103 (Mean): 
## 
##   Est.   S.E.   t val.      p
## ------ ------ -------- ------
##   0.03   0.04     0.83   0.41
## 
## Slope of mri_age_y when TissueFraction_GMWM = 0.7501615 (+ 1 SD): 
## 
##    Est.   S.E.   t val.      p
## ------- ------ -------- ------
##   -0.04   0.06    -0.74   0.46
#save Johnson-Neyman plot
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_ins_age_tf_int_JNplot.png")
## Saving 7 x 5 in image
#create interaction plot (make it pretty later)
#basic plot
Ins_age_tf_int_plot <- interact_plot(lag_ins_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM)
Ins_age_tf_int_plot

#linearity check
Ins_age_tf_int_plot_lin <- interact_plot(lag_ins_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM, modx.values="terciles", plot.points=TRUE, linearity.check=TRUE)
## Medians of each tercile of TissueFraction_GMWM are 0.539, 0.655, 0.747
Ins_age_tf_int_plot_lin
## `geom_smooth()` using formula 'y ~ x'

##Print model summaries Print formatted model summaries and tables using report and sjplot package

library(report)
#library(sjplot) #try for pretty summary tables

summary(lag_naa_age_TisFrac)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    858.0    876.8   -424.0    848.0      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9925 -0.5718 -0.0199  0.5594  3.9440 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1630   0.4037  
##  Residual             0.7184   0.8476  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)          11.22872    0.43459 317.97211  25.838  < 2e-16 ***
## mri_age_y             0.17233    0.02961 317.97750   5.819 1.44e-08 ***
## TissueFraction_GMWM   3.35162    0.51255 313.38178   6.539 2.51e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.703       
## TssFrc_GMWM -0.915  0.394
report(lag_naa_age_TisFrac)
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## We fitted a linear mixed model (estimated using ML and nloptwrap optimizer) to predict NAA_conc_molal with mri_age_y and TissueFraction_GMWM (formula: NAA_conc_molal ~ mri_age_y + TissueFraction_GMWM). The model included subj_id as random effect (formula: ~1 | subj_id). The model's total explanatory power is substantial (conditional R2 = 0.30) and the part related to the fixed effects alone (marginal R2) is of 0.15. The model's intercept, corresponding to mri_age_y = 0 and TissueFraction_GMWM = 0, is at 11.23 (95% CI [10.37, 12.08], t(313) = 25.84, p < .001). Within this model:
## 
##   - The effect of mri age y is statistically significant and positive (beta = 0.17, 95% CI [0.11, 0.23], t(313) = 5.82, p < .001; Std. beta = 0.32, 95% CI [0.21, 0.43])
##   - The effect of TissueFraction GMWM is statistically significant and positive (beta = 3.35, 95% CI [2.34, 4.36], t(313) = 6.54, p < .001; Std. beta = 0.35, 95% CI [0.25, 0.46])
## 
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using
report_table(lag_naa_age_TisFrac)
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## Parameter           | Coefficient |         95% CI | t(313) |      p | Effects |    Group | Std. Coef. | Std. Coef. 95% CI |    Fit
## -----------------------------------------------------------------------------------------------------------------------------------
## (Intercept)         |       11.23 | [10.37, 12.08] |  25.84 | < .001 |   fixed |          |       0.02 |     [-0.11, 0.14] |       
## mri age y           |        0.17 | [ 0.11,  0.23] |   5.82 | < .001 |   fixed |          |       0.32 |     [ 0.21, 0.43] |       
## TissueFraction GMWM |        3.35 | [ 2.34,  4.36] |   6.54 | < .001 |   fixed |          |       0.35 |     [ 0.25, 0.46] |       
##                     |        0.40 |                |        |        |  random |  subj_id |            |                   |       
##                     |        0.85 |                |        |        |  random | Residual |            |                   |       
##                     |             |                |        |        |         |          |            |                   |       
## AIC                 |             |                |        |        |         |          |            |                   | 857.96
## BIC                 |             |                |        |        |         |          |            |                   | 876.77
## R2 (conditional)    |             |                |        |        |         |          |            |                   |   0.30
## R2 (marginal)       |             |                |        |        |         |          |            |                   |   0.15
## Sigma               |             |                |        |        |         |          |            |                   |   0.85
summary(lag_cre_age_TisFrac)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    760.2    779.0   -375.1    750.2      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7493 -0.5056  0.0270  0.4419  4.8839 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08547  0.2923  
##  Residual             0.54120  0.7357  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           7.04677    0.36923 319.70413  19.085   <2e-16 ***
## mri_age_y             0.05661    0.02497 319.62390   2.267   0.0241 *  
## TissueFraction_GMWM   4.47115    0.43745 317.80861  10.221   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y
## mri_age_y   -0.701       
## TssFrc_GMWM -0.918  0.393
report(lag_cre_age_TisFrac)
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## We fitted a linear mixed model (estimated using ML and nloptwrap optimizer) to predict Cre_PCr_conc_molal with mri_age_y and TissueFraction_GMWM (formula: Cre_PCr_conc_molal ~ mri_age_y + TissueFraction_GMWM). The model included subj_id as random effect (formula: ~1 | subj_id). The model's total explanatory power is substantial (conditional R2 = 0.35) and the part related to the fixed effects alone (marginal R2) is of 0.25. The model's intercept, corresponding to mri_age_y = 0 and TissueFraction_GMWM = 0, is at 7.05 (95% CI [6.32, 7.77], t(315) = 19.08, p < .001). Within this model:
## 
##   - The effect of mri age y is statistically significant and positive (beta = 0.06, 95% CI [7.48e-03, 0.11], t(315) = 2.27, p = 0.024; Std. beta = 0.12, 95% CI [0.02, 0.22])
##   - The effect of TissueFraction GMWM is statistically significant and positive (beta = 4.47, 95% CI [3.61, 5.33], t(315) = 10.22, p < .001; Std. beta = 0.53, 95% CI [0.43, 0.63])
## 
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using
summary(lag_cho_age_TisFrac_int)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     23.3     46.0     -5.7     11.3      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5480 -0.6259 -0.0670  0.5067  4.4339 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01168  0.1081  
##  Residual             0.05179  0.2276  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                                Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                     2.45851    0.26231 312.15748   9.373   <2e-16
## mri_age_y                       0.03027    0.03614 303.65993   0.838    0.403
## TissueFraction_GMWM             0.01894    0.40778 309.31482   0.046    0.963
## mri_age_y:TissueFraction_GMWM  -0.11820    0.05899 301.29296  -2.004    0.046
##                                  
## (Intercept)                   ***
## mri_age_y                        
## TissueFraction_GMWM              
## mri_age_y:TissueFraction_GMWM *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y TF_GMW
## mri_age_y   -0.943              
## TssFrc_GMWM -0.981  0.948       
## m__:TF_GMWM  0.896 -0.976 -0.942
report(lag_cho_age_TisFrac_int)
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## 'interpret_d()' is now deprecated. Please use 'interpret_cohens_d()'.
## We fitted a linear mixed model (estimated using ML and nloptwrap optimizer) to predict Cho_GPC_PCh_conc_molal with mri_age_y and TissueFraction_GMWM (formula: Cho_GPC_PCh_conc_molal ~ mri_age_y * TissueFraction_GMWM). The model included subj_id as random effect (formula: ~1 | subj_id). The model's total explanatory power is substantial (conditional R2 = 0.28) and the part related to the fixed effects alone (marginal R2) is of 0.12. The model's intercept, corresponding to mri_age_y = 0 and TissueFraction_GMWM = 0, is at 2.46 (95% CI [1.94, 2.97], t(314) = 9.37, p < .001). Within this model:
## 
##   - The effect of mri age y is statistically non-significant and positive (beta = 0.03, 95% CI [-0.04, 0.10], t(314) = 0.84, p = 0.403; Std. beta = -0.32, 95% CI [-0.44, -0.21])
##   - The effect of TissueFraction GMWM is statistically non-significant and positive (beta = 0.02, 95% CI [-0.78, 0.82], t(314) = 0.05, p = 0.963; Std. beta = -0.27, 95% CI [-0.38, -0.16])
##   - The interaction effect of TissueFraction GMWM on mri age y is statistically significant and negative (beta = -0.12, 95% CI [-0.23, -2.13e-03], t(314) = -2.00, p = 0.046; Std. beta = -0.09, 95% CI [-0.18, -1.62e-03])
## 
## Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using

##Plot results of models including tissue fraction Estimate metabolite residuals using remef Include individual slopes (estimated via linear models by subject) using geom_smooth() Include overall fit line using geom_abline

#LAG - NAA

#pull fixed effects for overall model
fixef.lag_naa_TisFrac<-fixef(lag_naa_age_TisFrac)

#subset dataframe with only subjects included in NAA analysis
lag_dat_plot<- lag_dat %>% 
  filter(! is.na(NAA_conc_molal))

#calculate partial effect of age on tNAA (remove effect of tissue fraction)
lag_naa_partial <- remef(lag_naa_age_TisFrac, fix = "TissueFraction_GMWM") #does NOT remove random effects variance

#create plot to show NAA ~ age while controlling for Tissue Fraction
plot_lag_naa_age_tf<-ggplot(lag_dat_plot, aes(x=mri_age_y, y=lag_naa_partial, na.exclude = TRUE, color=as.factor(female))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, aes(x=mri_age_y, y=lag_naa_partial, group=as.factor(subj_id)), size = .2, show.legend = FALSE) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA | Tissue Fraction") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_naa_TisFrac[[1]], slope = fixef.lag_naa_TisFrac[[2]])) + 
  theme_classic()
plot_lag_naa_age_tf
## `geom_smooth()` using formula 'y ~ x'

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_naa_age_partialef.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
#LAG - Cre

#pull fixed effects for overall model
fixef.lag_cre_age_TisFrac<-fixef(lag_cre_age_TisFrac)

#subset dataframe with only subjects included in NAA analysis
lag_dat_plot_Cre<- lag_dat %>% 
  filter(! is.na(Cre_PCr_conc_molal))

#calculate partial effect of age on Cre (remove effect of tissue fraction)
lag_cre_partial <- remef(lag_cre_age_TisFrac, fix = "TissueFraction_GMWM") #does NOT remove random effects variance

plot_lag_cre_age_tf<-ggplot(lag_dat_plot_Cre, aes(x=mri_age_y, y=lag_cre_partial, na.exclude = TRUE, color=as.factor(female))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, aes(x=mri_age_y, y=lag_cre_partial, group=as.factor(subj_id)), size = .2, show.legend = FALSE) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tCre | Tissue Fraction") +
  scale_color_viridis(name = "Sex", labels = c("Male", "Female"), discrete = TRUE, option = "H", begin = .2, end = .8) +
  geom_abline(aes(intercept = fixef.lag_cre_age_TisFrac[[1]], slope = fixef.lag_cre_age_TisFrac[[2]])) + 
  theme_classic()
plot_lag_cre_age_tf
## `geom_smooth()` using formula 'y ~ x'

  ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_cre_age_partialef.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
#LAG - Cho
Cho_age_tf_int_plot <- interact_plot(lag_cho_age_TisFrac_int, pred = mri_age_y, modx = TissueFraction_GMWM, data=lag_dat, plot.points=TRUE, legend.main="Tissue Fraction GM/(GM+WM)", modx.labels=c(".535 (-1 SD)", ".642 (Mean)", ".750 (+1 SD)")) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tCho (molal)") + 
  scale_color_viridis(option="mako", name="Tissue Fraction GM/(GM+WM)", begin= .32, end=1) +
  theme_classic()
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
Cho_age_tf_int_plot
## Warning: Removed 1 rows containing missing values (geom_point).

  ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_cho_age_tf_int.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).